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A VORTEX- LATTICE METHOD FOR THE MEAN CAMBER SHAPES 
OF TRIMMED NONCOPLANAR PLANFORMS 
WITH MINIMUM VORTEX DRAG 

John E. Lamar 
Langley Research Center 

SUMMARY 

A new subsonic method has been developed by which the mean camber surface can 
be determined for trimmed noncoplanar planforms with minimum vortex drag. This 
method uses a vortex lattice and overcomes previous difficulties with chord loading 
specification. This method uses a Trefftz plane analysis to determine the optimum span 
loading for minimum drag, then solves for the mean camber surface of the wing, which 
will provide the required loading. Pitching- moment or root-bending-moment constraints 
can be employed as well at the design lift coefficient. 

Sensitivity studies of vortex-lattice arrangement have been made with this method 
and are presented. Comparisons with other theories show generally good agreement. 

The versatility of the method is demonstrated by applying it to (1) isolated wings, (2) wing- 
canard configurations, (3) a tandem wing, and (4) a wing-winglet configuration. 

INTRODUCTION 

Configuration design for subsonic transports usually begins with the wing, after 
which the body and its effects are taken into account, and then the tails are sized and 
located by taking into account stability and control requirements. With the advent of 
highly maneuverable aircraft having closely coupled lifting surfaces, there has been an 
increased interest in changing the design order so that multiple surfaces could be 
designed together to yield a trimmed configuration with minimum induced drag at some 
specified lift coefficient. Such a combined design approach requires that the mutual inter- 
ference of the lifting surfaces be considered initially. 

Single planform design methods are available to optimize the mean camber surface, 
better called the local elevation surface, for wings flying at subsonic speeds (for example, 
ref. 1) and at supersonic speeds (for example, refs. 2 and 3). The design method pre- 
sented in reference 1 was developed from an established analysis method (Multhopp type), 
also presented in reference 1, by using the same mathematical model, but the design 


method solves for the local mean slope rather than the lifting pressures. In the usual 
implementation of reference 1, the design lifting pressures are taken to be linear chord- 
wise, but must be represented in this solution by a sine series which oscillates about 
them. An example presented herein demonstrates that corresponding oscillations may 
appear in pressure distributions measured on wings which have been designed by the 
method of reference 1. The method developed herein overcomes this oscillatory lifting 
pressure behavior by specifying linear chord loadings at the outset. 

The development approach used in the two-planform design problem will be simi- 
lar to that used for a single planform. The analytic method employed, selected because 
of its geometric versatility, is the noncoplanar two-planform vortex-lattice method of 
reference 4. 

The design procedure is essentially an optimization or extremization problem. 
Subsonic methods (for example, see refs. 5 and 6) are available for determining the span 
load distributions on bent lifting lines in the Trefftz plane, but they do not describe the 
necessary local elevation surface. This is one of the objectives of the present method 
which will utilize the Lagrange multiplier technique (also employed in refs. 2 and 3). The 
method of reference 4 is used to provide the needed geometrical relationships between the 
circulation and induced normal flow for complex planforms, as well as to compute the lift, 
drag, and pitching moment. 

This paper also presents the results of precision studies and comparisons with 
other methods and data. Several examples of solutions for configurations of recent inter- 
est are also presented. The FORTRAN computer program written to perform the compu- 
tation is described (appendix A), along with details of the program input data (appendix B) 
and output data (appendix C). Listings and typical running times of example configura- 
tions are given (appendix D), and a FORTRAN program listing is provided (appendix E). 
Appendix F provides details concerning the changes needed to substitute a root-bending- 
moment constraint for the basic constraint on configuration pitching- moment balance. 

SYMBOLS 

The geometric description of planforms is based on the body-axis system. (See 
fig. 1 for positive directions.) For computational purposes the planform is replaced by a 
vortex lattice which is in a wind-axis system. Both the body axes and the wind axes have 
their origins in the planform plane of symmetry. (See sketch (a) for details.) The axis 
system of a particular horseshoe vortex is wind oriented and referred to the origin of that 
horseshoe vortex (fig. 1). For the purpose of the computer program, the length dimension 
is arbitrary for a given case; angles associated with the planform are always in degrees. 
(The variable names used for input data in the computer program are described in 
appendix B.) 
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element of influence function matrix A, 
tains induced normal flow at Zth point 
strength; total number of elements is 


■ ^v,Z,n 


4 tt 


, which con- 


due to nth horseshoe vortex of unit 

N V N 
2 2 


AR aspect ratio 


a 


t 

I b 



fractional chord location where chord load changes from constant value to 
linearly varying value toward zero at trailing edge 

coefficients in spanwise scaling polynomial 


wing span 


root -bending- moment coefficient about X-axis, 


Root bending moment 


drag coefficient, 


Drag 

^oo^ref 


Cd q drag coefficient at = 0 


lift coefficient. 


Lift 

^oo^ref 


^m 


pitching- moment coefficient about Y-axis, 


Pitching moment 
^^oc^ref^ref 


'N 


normal-force coefficient. 


Normal force 
^oo^ref 


ACp lifting pressure coefficient 

c chord 


c^ section lift coefficient 

Cref reference chord 
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influence function which geometrically relates induced effect of nth 
horseshoe vortex to quantity which is proportional to induced down- 
wash or sidewash at slope point I (see sketch (a) and also eqs. (5) 
and (6)) 

sum of influence function ^ or ^ at slope point I on 

planform caused by two symmetrically located horseshoe vortices, 
left wing panel vortex denoted by n and right wing panel vortex 
denoted by N + 1 - n (see fig. 1) 

function to be extremized (see eq. (19)) 

” — I 

= Nca + 0.75 ! (brackets indicate "take the greatest integer") 


maximum number of spanwise scaling terms (see eqs. (25) to (27)) 


pitching moment about coordinate origin 

free -stream Mach number 

number of span stations where pressure modes are defined as used in 
reference 1 

maximum number of elemental panels on both sides of configuration; maxi- 
mum number of chordal control points at each of m span stations as used 
in reference 1 

number of elemental panels from leading to trailing edge in chordwise row 

total number of (chordwise) rows in spanwise direction of elemental panels 
on configuration semispan 

free -stream dynamic pressure 


reference area 
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horseshoe vortex semiwidth in plane of horseshoe (see fig. 2) 


i ^ 

X,Y,Z 
I X,Y,Z 

i x,iz 

I 

x,y,z 

x’ 

! 

I x,y,z 
Ax 

^ ^c/4 

^3c/4 

y*,z* 


Zc 

z/c 


(8z/9x)^ 

a 




n 


free -stream velocity 

axis system of given horseshoe vortex (see fig. 1) 
body-axis system for planform (see fig. 1) 
wind-axis system for planform (see sketch (a)) 
distance along X-, Y-, and Z-axis, respectively 

= x//3 

distance along X-, Y-, and Z-axis, respectively 

incremental movement of X-Y coordinate origin in streamwise direction 

midspan x-location of quarter-chord of elemental panel 

midspan x-location of three-quarter-chord of elemental panel 

y and z distances from image vortices located on right half of plane of 
symmetry, as viewed from behind, to points on left panel 

canard height with respect to wing plane, positive down 

local elevation normalized by local chord, referenced to local trailing-edge 
height, positive down 

ith elemental local slope in vector {9z/8x^ of N/2 elements (see eq. (1)) 
angle of attack, deg 

Prandtl-Glauert correction factor to account for effect of compressibility in 
subsonic flow, yl - 

vortex strength of nth element in vector of N/2 elements 
independent variable in extremization process 
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e 


incidence angle, positive leading edge up, deg 


rj nondimensional spanwise coordinates, 

b/2 

Vi nondimensional spanwise coordinate based on local planform semispan 

A planform leading-edge sweep angle in X-Y plane, deg 

X Lagrange multiplier (see eq. (19)) 

I distance along local chord normalized by local chord 

^ fractional chordwise location of point where mean camber height is to be 

computed (see eq. (28)) 

dihedral angle from trailing vortex to point on left panel being influenced; 
a measured from left panel, o’ measured from right panel 

0 constraint function (see eqs. (20) and (21)); also horseshoe vortex dihedral 

angle in Y-Z plane on left wing panel, deg 

0' horseshoe vortex dihedral angle on right wing panel, cp' = -(p, deg 

0 quarter-chord sweep angle of elemental panel; because of small angle 

assumption, also used as sweep angle of spanwise horseshoe vortex fila- 
ment in X-Y plane, deg 

=tan-l(^j 

Subscripts: 
c canard 

d design 

i,j,k indices to vary over the range indicated 

le leading edge 
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I 


i associated with slope point and horseshoe vortex, respectively, ranging from 

1 to N/2 

L left trailing leg 

R right trailing leg 

r root-chord location 

V vortex 

1 w wing 

' Matrix notation: 


O 


column vector 


' [ ] 
L 


square matrix 


X 



u 



Flow angle of attack determined 
at each slope point 

Typical spanwise 
vortex filament 




Vortex -lattice trailing filaments 



Sketch (a) 


THEORETICAL DEVELOPMENT 

This section presents the application of vortex-lattice methodology to the mean- 
camber-surface design of two lifting planforms which may be separated vertically and 
have dihedral. For a given planform, local vertical displacements of the surfaces with 
respect to their chord lines in the wing axis (see sketch (a)) are assumed to be negligi- 
ble; however, vertical displacements of the solution surfaces due to planform separation 
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or dihedral are included. The wakes of these bent lifting planforms are assumed to lie in 
their respective extended bent chord planes with no roll up. For a two-planform configu- 
ration the resulting local elevation surface solutions are those for which both the vortex 
drag is minimized at the design lift coefficient and the pitching moment is constrained to 
be zero about the origin. For an isolated planform no pitching-moment constraint is 
imposed. Thus, the solution is the local elevation surface yielding the minimum vortex 
drag at the design lift coefficient. Lagrange multipliers together with suitable interpo- 
lating and integrating procedures are used to obtain the solutions. The details of the 
solution are given in the following five subsections. 

Relationship Between Local Slope and Circulation 

From reference 4, the distributed circulation over a lifting system is related to the 
local slope by 



where the matrix [aJ is the aerodynamic influence coefficient matrix based on the paneling 
technique described in reference 4. This matrix has elements of 


‘ ><t>) tan 0^ 


( 2 ) 


which, because of the assumed spanwise symmetry of loading, leads to 


panel 


" panel 


(3) 


and 




^v,Z,N+l-n^^ >y>z,s,0 panel 


(4) 
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where 


(y tan - x') cos <t> 

F^(x ,y,z,s,i^ ,(l>) = j 1 r ^ 

(x') + (y sin <p) + cos2 0(^y2 tan2 \p' + z2 sec2 \p' - 2yx' tan \p') - 2z cos 0 sin 0(y + x' tan 


(x' + s cos (b tan \p') cos 0 tan ip' + {y + s cos <b) cos <^ + (z + s sin <p) sin 0 


-jl/2 

(x’ + s cos 0 tan 0')^ + (y + s cos 0)^ + (z + s sin 0)^ 


(x* - s cos 0 tan 0') cos 0 tan 0' + (y - s cos 0) cos 0 + (z - s sin 0) sin 0I 


t 2 2 2 

(x' - s cos 0 tan 0 ) + (y - s cos 0) + (z - s sin 0) 


nl/2 


y - s cos 0 


(y - s cos 0)^ + (z - s sin 0)^ 


- • 


x’ - s cos 0 tan 0' 


(x’ - s cos 0 tan 0’)^ + (y - s cos 0)^ + (z - s sin 0)^ 


y + s cos 0 


2 2 
(y + s cos 0) + (z + sin 0) 


x’ + s cos 0 tan 0' 


2 2 2 
(x’ + s cos 0 tan 0 ) + (y + s cos 0) + (z + s sin 0) 


ll/2f 


( 5 ) 


and 


Fy(x’,y,z,s,i//',(i>) = 


x’ sin 0 - z cos 0 tan 0* 


(x')^ + (y sin 0)^ + cos2 0(y2 tan2 0' + z2 sec2 0' - 2yx' tan 0’) - 2z cos 0 sin 0(y + x’ tan 0') 


J(x' + s cos 0 tan 0’) cos 0 tan 0’ + (y + s cos 0) cos 0 + (z + s sin 0) sin 0 


[ 2 2 2I 

(x’ + s cos 0 tan 0 ) + (y + s cos 0) + (z + s sin 0) 


1/2 


(x' - s cos 0 tan 0') cos 0 tan 0' + (y - s cos 0) cos 0 + (z - s sin 0) sin 0l 


2 2 2I 

(x’ - s cos 0 tan 0’) + (y - s cos 0) + {z - s sin 0) 


,1/2 


z - s sin 0 


x’ - s cos 0 tan 0’ 


2 2 
(y - s cos 0) + (z - s sin 0) 


(x’ - s cos 0 tan 0’)"^ + (y - s cos 0)"^ + (z - s sin 0) 


2I 


1/2 


z + s sin 0 


1 - ■ 


X* + s cos 0 tan 0 ^ 


(y + s cos 0) + (z + s sin 0) ^ g ^ _j_ g ^^2 ^ + s sin 0)^1 


1/2 


( 6 ) 


with I signifying the particular slope point and n the particular horseshoe vortex 
influencing the slope point. 
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Circulation Specification 

Once the surface slope matrix (dz/ 8x} is known, chordwise integration can be per- 
formed to determine the local elevation surface z/c, which contains the effects of camber, 
twist, and angle of attack. The major problem to be solved is determining the necessary 
circulation matrix {^r/U} to employ in equation (1). The problem is simplified somewhat 
by having the chordwise shape of the bound circulation remain unchanged across each span, 
although the chordwise shape may vary from one planform to another. (This simplification 
can easily be removed without any new analysis and would require only a small programing 
change.) The chordwise loadings allowable in the program range from rectangular to right 
triangular toward the leading edge and were selected because they are of known utility. An 
example is given in sketch (b). Two different techniques are utilized to arrive at the span- 
wise scaling of the chordwise shapes. The particular technique to be employed depends on 
whether the configuration has dihedral. 


Chordwise (typical) Assumed spanwise modes 

(configurations without dihedral) 



Sketch (b) 
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For a configuration having dihedral, the spanwise scaling must be determined dis- 
cretely because no finite polynomial representation of the scaling is known with certainty, 
even for an isolated wing. However, for configurations with no dihedral, the spanwise 
scaling can be written as a polynomial for each planform, 


(see sketch (b)) with a maximum of three coefficients per planform being determined as 
part of the solution. It is possible to write this polynomial as a solution because the iso- 
lated wing solution is known to be of the elliptical form \jl - and the presence of the 

other planform is assumed to generate a loading disturbance which can be represented by 
the other two terms in addition to adjusting a^. Once the scaling is known from either 
technique, then {t/U) is readily obtained by multiplication. 


Lift, Pitching-Moment, and Drag Contributions 

The contributions to Cj^ and to Cm, respectively, from the jth chordwise row of 
horseshoe vortices are 


and 


Cr . = 


q=oSref 


l_ ^ cos 
'ref ^oo^ref 


( 7 ) 


My ,3 _ 4q^s cos 0 


*loo^ref^ref ^oo^ref^ref 


‘hh 


( 8 ) 


where 


U/i 




1-^i 


1 - a 






li^a 


(,>a 


(9 a) 


t _ i - 0.75 


(9b) 
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It should be observed that no contribution from the drag forces is included in equation (8). 

I 

Even though . and Cj^ i actually occur on the wing at the jth spanwise loca- 

tion, they can be utilized in a Trefftz plane solution if the chordwise summations are per- 
formed. This utilization is possible herein because the trailing wake is assumed not to 
roll up, and the general configuration has specifiable chord loading shapes. Summing the 
chordwise loadings at this point allows the solution of the spanwise scaling to be per- 
formed on a bent lifting line located in the Trefftz plane, which is, of course, ideally 
suited for the vortex drag computation. In addition, the summation reduces the number 
of unknowns from the product of N(> and Ng to only Ng. Hence, a larger value of Ng 
can be used in the Trefftz plane, which should yield improved accuracy in the spanwise 
scaling factors without affecting the number of horseshoe vortices on the wing. Then, 
when the circulations are needed on the wing for use in equation (1), the well-defined var- 
iations of the spanwise scaling factors are interpolated to the original spanwise positions | 
of the wing vortex lattice which is used to generate [a]. The procedure is implemented 

as follows: 

The summation in the lift expression (eq. (7)) can be written as 



where I is the last i value which satisfies = a; that is, 


I = 


Nca 


+ 0.75 


where the brackets indicate "take the greatest integer." Hence, 


(12) 



= I + 


(Nc + 0.75)(Nc - I 
N„(l - a) 


N, 


^ 1 - y i 

Nc(l - a) 


( 13 ) 


Similarly, the summation in the pitch expression (eq. (8)) can be written as 





L 


0.75cj (Nc + 0.75)(Nc - ij 

^Tjl ^ VI -a) 


I 


1 


/ i 

i=l 


N„(l - a) 


(^le)j 


1.5C; 




Nr 


Nc 

, .) 

i ik+1 


Nc 

N„ (1 - a) i=l+i 


"j 


( 14 ) 


The contribution to the vortex drag coefficient at the ith chordwise row due to the 
jth chordwise row is obtained by using only half the trailing vortex induced normal wash 
from the Trefftz plane. The result is 



In the ± sign, plus indicates that the trailing vortex filament is to the left of the influenced 
point; minus, to the right. 

In using equations (7), (8), and (15), a new vortex system is set up in the Trefftz 
plane in which the bent chord plane is represented by a system of uniformly spaced trail- 
ing vortices (the quantity 2s in fig. 2). This uniformity of vortex spacing leads to a 
simplification in the equations and can be thought of as a discretization of the ideas of 
Munk (ref. 7) and Milne -Thomson (ref. 8) for a bound vortex of constant strength. 
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Spanwise Scaling Determination 

To determine the spanwise scaling with either technique requires the combination 

of the contributions from each spanwise position for configurations with dihedral or the 

mode shape contributions for configurations without dihedral. These contributions must 

be employed in the appropriate total Cl and Cm constraint equations as well as in the 

Cr. extremization operation. The details of the solution for configurations with dihe- 
u,v 

dral are as follows: 


and 


j=l 

Ns 

Cm = 2 ^ 5frn,i 
j=l 

Ns Ns 
i=l j=l 


( 16 ) 


(17) 


(18) 


where the terms are the spanwise scaling factors and the independent variables in 
the solution. 

The problem is formalized in the Lagrange extremization method by forming the 
fmiction to be extremized 


2 

^ " ^D,v Jj ^i'^i 
i=l 

with the two constraint equations 
Ns 

^1=^ Ij ®k^L,k ■ ^L,d ^ ® 
k=l 

^2 ~ ^ ®k^m,k ”0 = 0 
k=l 


(19) 


( 20 ) 


( 21 ) 
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where Xj and X 2 are the Lagrange multipliers. In order to extremize the function G, 
it is necessary to find a solution to the set of linear equations resulting from 

(f = 1, 2, . . Ns) (22) 

and 

— = 0 (i = 1, 2) (23) 



where equation (23) is just a restatement of equations (20) and (21), The Ng equations 
represented in equation (22) are explicitly 

Ng 

^ (^D,i,k + ^D,k,i)^k + ^L.i^l ^m,i^2 = ® (i = 1, 2, . . Ng) (24) 

k=l ^ 

Equations (24), (20), and (21) provide Ng + 2 relations having as the Ng + 2 unknowns 
the Ng values of Xj^, and X 2 . 

The matrix to be solved for configurations with dihedral can be as large as 102 
square, and it is possible for this matrix to become ill conditioned if the trailing vortex 
filaments from the two planforms coincide. If this coincidence occurs, an alternative 
matrix inversion routine, based on least squaring, is utilized. 

It is difficult to assess the accuracy of the calculated values of 5jj because mini- 
mum vortex drag y) solutions are not generally known, even for isolated wings having 

dihedral. As a numerical check, the ratio of the normal induced velocity to the cosine of 
the local dihedral angle is computed. According to Munk (ref, 7), this ratio should be con- 
stant across the configuration span for minimum vortex drag. Hence, the uniformity of 
this ratio is an indication of the accuracy of those solutions for which only the lift con- 
straint is operative. If both the lift and moment constraints are operative, then the vor- 
tex drag will be the minimum obtainable for the problem posed but not necessarily an 
absolute minimum. Under the pitching-moment-constraint conditions, this numerical 
check is meaningless and should be ignored. 

It should be noted that equation (21) could be changed from a pitching- moment con- 
straint to one which involved the root bending moment. In fact, this has been done in one 
of the examples discussed in the text. Details for implementing this constraint are given 
in appendix F, 
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where K = 6 and Cj^ and Cm are the and contributions associated 

with the kth term in the polynomials 






(Note that k = 1, 2, and 3 are assigned to the first planform and 4, 5, and 6 to the sec- 
ond.) These contributions are computed by first assuming a unit value of scaling with 
each term in the polynomial, then multiplying each resulting spanwise scaling distribution 
by the Cj^ ^ and Cm^j terms of equations (7) and (8), and finally summing spanwise 

over all the chordwise rows associated with each set of k values (or planform). The 
vortex drag coefficient associated with the ith and kth combination of spanwise scaling 
distributions is computed similarly. The 6jj terms are equivalent to the 

unknown coefficients in the polynomial and are the independent variables in the solution. 

The extremization of equation (27), with the same Cj^ and Cm constraints as 
before, produces K + 2 relations with the K values of 6jj, Xj, and X2 as the 
unknowns. Obviously, this matrix, no larger than 8 square, is much smaller and hence 
faster to invert than that utilized for configurations having dihedral. 
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Determination of Local Elevation Curves 


With 5jj known, then {t/U}, Cj^, Cm, and Cq ^ can be determined. The 

results for {t/u} are interpolated to the original spanwise positions of the paneling 
which is used in equation (1) and in the following equation to find the local elevation 
curves. The equation for the local elevation above the computational plane at a 
particular point (4',y) is 

t* _ 

I (^’,y) = f |la,y) ( 28 ) 

^ o'! 8x 


Further discussion is given regarding this integration in the section "Precision," but it 
should be noted that cubic splines are utilized to interpolate the local surface slopes 
between slope points as well as to integrate the resulting distribution. 

Three additional aspects of the present method should be noted: (1) The local slope 
and elevation results obtained are linearly dependent on Cj_^; hence, they can be used to 
obtain design information at other than the original design Cj^ by multiplying these 
results by the ratio of the new value to the old value of Cj_^. (2) For an isolated plan- 
form with zero dihedral, the three assumed spanwise distributions are self-reducing; 
that is, the Lagrange multipliers of the second and third distributions become zero, leav- 
ing only the first (the elliptic form) to give the correct minimum vortex drag. Thus, only 
the elliptic spanwise distribution is imposed for mean-camber-surface solutions of iso- 
lated planforms without dihedral. (3) As a result of the relationship between F/U and 
the lift on an elemental panel, F/U is related to the assumed constant value of ACp 
over the panel by 


F_ ^^P c 
U" 2 Ne 


(29) 


for a viniform chordwise distribution of elemental panels. If a nonuniform distribution is 
used, then equation (29) as well as the computer implementation must be modified. 


RESULTS AND DISCUSSION 


General 

Before the design method just outlined is employed, it is necessary to examine the 
sensitivity of its results to vortex-lattice arrangement. It is also important to compare 
results obtained with this method with those available in the literature. Unfortunately, 
the available solutions, whether exact or numerical, may not be for configurations which 
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will exercise the constraint or extreminization capabilities of the present method. In 
fact, the available exact solutions are for configurations which are either two-dimensional 
sections or isolated three-dimensional wings with a nonelliptic span loading. The solu- 
tions for such configurations require program modifications to the span loading and 
involve no optimization. The numerical solution used for comparison is for an isolated 
planform without dihedral. 

Two-dimensional comparisons are used to determine suitable chordwise locations 
and the number N^. of horseshoe vortices. The effect that different extrapolations of 
the chordwise representation of 9z/ 9x ahead of the first and behind the last slope points 
have on the local elevation curve has also been investigated. In addition, the sensitivity 
of the local elevation solutions to the number Ng and location of chordwise rows of 
horseshoe vortices was investigated for an isolated planform. 

Following the section "Precision," an application of the present method to a wing- 
canard configuration is given for various vertical separations and moment trim points, as 
well as a comparison of the local incidence distributions, vortex drag values, and span 
loadings. Calculated results for a tandem wing and for a wing-winglet combination are 
also presented. 


Precision 

Results of the present method are presented in figures 3, 4, 5, 6, and 7 in terms of 
local midsurface slopes and elevations along the chord. The local elevation results can 
be thought of as including the effects of incidence, twist, and mean camber. In these fig- 
ures results of the present method are compared with those obtained from other methods. 
Where appropriate, these comparisons are made at a number of spanwise locations. 

Two-dimensional .- Figure 3 presents, for the three particular chordwise lifting 
pressure distributions shown in the inset sketches, the local slope and elevation varia- 
tions along the chord obtained from the present method and from two-dimensional theory 
(ref. 9) at ^ = 1.0 and Moo = 0. The predictions of the present method were obtained 
by utilizing an aspect-ratio- 50 rectangular wing with a rectangular span loading. In order 
to avoid the tip effects, only the results near the plane of symmetry are presented and they 
are for the different chordwise patterns and extrapolations denoted. It should be noted that 
the diamond symbol does not appear in those parts of figure 3 which present the local slope 
since the results are coincident with the results denoted by the square symbol. Only in the 
parts of the figure presenting local elevation do the two symbols not coincide, which is a 
result of the differences evolving from methods of extrapolation. A discussion of the 
extrapolation methods will be presented later in this section. 

The local elevation curves are identified by values of a of 0.2, 0.6, and 1.0 which 
denote the fractional chordwise locations where the net pressure changes from a constant 
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value to a linearly varying value toward zero. (The symbol a employed in the text and 
figures is the same as the variables XCFW or XCFT used in appendix B.) 

The following observations can be made from figure 3; 

(1) For all patterns studied and for all values of a, the present method predicts 
both the amplitude and trend of the local slope analytic curves with reasonable accuracy. 

(2) Of the uniform chordwise vortex patterns, the one with = 20 is superior to 
that with Np = 10. This result can be attributed to two causes: (a) increasing N(. 
from 10 to 20 provides more definition to the approximate curve, especially near the 
chordwise edges, where the analytic result may have a steep gradient, and (b) the extrap- 
olations to the chordal edges, which must be employed with the approximate curve for 
integration purposes, are more accurate as a result of the smaller distance over which 
they must be applied. Because the chordwise integration of local slope occurs from the 
trailing edge forward, any errors in matching the analytic local slopes at or ahead of the 
trailing edge will be seen forward of that chordwise location and will accumulate. 

The incidence angles have been extracted from the local elevation curves for com- 
parison with the exact solution. The following table summarizes the incidence angles 
obtained with N^. = 10 and 20, as well as those for Ng > 20 from a modified version of 
the program. The results from > 20 are provided so that the solution convergence 
and its rate can be examined. The table clearly shows that the results of the present 
method are more positive than, but tend toward, the exact ones with increasing values of 
Nc but at a slower rate as N(. increases. Though not shown herein, it was observed 
that the results obtained for N(> = 40 also gave the best agreement with the local eleva- 
tion curve, especially for a = 1.0. Also, with Ng held constant, doubling the value of 
Nc provides a reduction in percent error of less than a factor of two while approximately 
quadrupling the computer time. Hence, N^. = 20 is the largest number that will be 
employed; however, the best chordwise spacing of these boimd vortices is still to be 
determined. 


a 

Incidence angle, deg, from — 

Present method (uniform spacing) 

Exact 
solution 
(ref. 9) 

Nc = 10 

Nc = 20 

o 

CO 

tl 

Nc=40 

0.2 

5.3359 

4.9097 

4.7504 

4.6650 

4.1752 

.6 

3.5421 

3.2109 


3.0167 

2.6052 

1.0 

1.3863 

.8594 


.5386 

0 
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(3) Because the = 20 uniform solution has slope points nearer the chordal 
edges and yielded better local elevation curves than the solution for N(. = 10, it was 
anticipated that by arranging the locations of the elemental panels nearer the chordwise 
edges, as in a cosine manner, there could be improvements in the = 10 solution. 

The results of this change are as follows: (a) Better agreement with the analytic local 
slope curves near the leading edge is produced for all values of a and near the trailing 
edge for a = 1.0. (b) There is poorer agreement from 0.1 to 0.7 chord where the new 
local slopes are less than those of the uniform spacing and the analytic curve. The error 
accumulates to a larger overestimation of local elevation from 0.1 to 0.5 chord than for 
the other patterns, (c) As a result of the better local slope prediction near the leading 
edge, the local elevation predictions at the leading edge are better with the cosine spacing 
than for the uniform spacing solution with the same number of divisions. 

The incidence angles have also been computed for the solutions just discussed along 
with those for Nc = 20 which employed a cosine spacing. These results are summarized 
in the following table: 


a 

Incidence angle, deg, from - 

Present method 
(cosine spacing) 

Exact 
solution 
(ref. 9) 

Nc = 10 

Nc = 20 

0.2 

4.9836 

4.6650 

4.1752 

.6 

3.3594 

3.0567 

2.6052 

1.0 

1.1859 

.6474 

0 


This table shows the improvement in incidence angle prediction with increasing values of 
which were obtained for a cosine spacing. If these results are compared with those 
from the previous table, it can be seen that the percent error decreases with a change 
from a uniform spacing to a cosine spacing at a given value of N^,. However, the local 
elevation solutions with the cosine spacing are generally poorer when compared with the 
exact solutions over the midchord range than those with the uniform spacing. Since the 
local elevation surface is the primary purpose of the computation, the uniform spacing is 
utilized in the following numerical studies and applications and is also employed in the 
program. 

(4) The solution for N^. = 20 uniformly spaced horseshoe vortices gives the best 
overall results. Furthermore, the effect of changing the method of extrapolation of the 
local slope curve ahead of the first and behind the last slope points is not significant. 

(See sketch (c) for a comparison of the two extrapolation techniques near the leading edge.) 
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Consequently, the constant extrapolation method is employed because of the projected 
computer resource savings. 

— • — dz/dx solution and fairing 

Spline extrapolation 

Constant extrapolation 


dz/dx 


0 .1 .2 .3 

Sketch (c) 

Number of rows along semispan ^Ng^- The number of chordwise rows Ng needed 

on a semispan is studied in three dimensions by using a low-aspect-ratio (2.50) trapezoidal 
wing as an example planform. The local slopes and elevations along the chord are pre- 
sented in figure 4. The fixed parameters are rectangular chord loading (a = 1.0), elliptic 
span loading, ^ = 0.35, Moo = 0.40, and Nq = 20. The value of Ns is taken to be 
either 10 or 20, with both uniform and cosine distributions employed over the semispan. 
From the data of figure 4, little sensitivity is noted across the semispan to either the num- 
ber or the distribution of chordwise rows used in the solutions. The one exception is at 

— — = 0.95 with Ng = 20 and a cosine distribution, where a local elevation surface with 
b/2 

a larger reversal in the incidence near the tip occurs. The incidence reversal is so great 
that it is suspect; better results could be obtained from a smooth fairing of the inboard 
results to the tip. 
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An effect of changing Ng is that there is little benefit to be gained by using large 
values of Ng, except for the expanded number of local elevation curves tabulated in the 
computer output. This lack of benefit is associated with prescribing the span and chord 
loadings in advance in the design problem for the wing without dihedral, whereas they 
must be determined locally in the analysis problem. Hence, less sensitivity in the results 
is noted. For the example wing, a value of Ng = 10 was found to be sufficient; however, 
Ng values of this magnitude may not be large enough for other planforms and Mach num- 
bers. A second effect is that a cosine distribution across the span of chordwise rows 
does not improve the solutions and can in fact lead to poorer ones because of an unrea- 
sonable incidence distribution for spanwise locations too near the tip. Hence, a uniform 
distribution of chordwise rows is recommended and utilized herein. 

Precision of the solution for wings with dihedral .- It is useful to consider whether 
the type of studies conducted for wings with no dihedral needs to be repeated for wings 
with dihedral. Since the difference in technique is limited to the procedure for determin- 
ing the spanwise scaling, the results of the Nc study should be valid for both techniques. 
Concerning the Ng study, the spanwise scaling differences are restricted to the optimi- 
zation part of the program where repaneling occurs. In either solution the set of answers 
is evaluated or interpolated to the original paneling scheme for the computation of ^z/8x^. 
Hence, it is only necessary to determine whether there are enough discrete spanwise scal- 
ing values to obtain a good approximation to the functional form of the solution. Thus, the 
technique usually used only when dihedral is present was applied to the isolated flat wing 
of figure 4, and the interpolated span loading results for both techniques are presented in 
figure 5. The agreement is seen to be generally good except in the outer 10 percent of 
semispan. In that region the functional form has the largest variation and is more diffi- 
cult to represent discretely. However, the discrete solution di^ yield a constant value 
of normal velocity across the span, which is the proper result. The c^c interpolated 
results shown in figure 5 for wings without dihedral are a part of the original elliptical 
curve. In addition, the difference in between the two techniques is 0.0008. Com- 

paring Cq ^ with Cj^^y^AR shows that the difference due to the technique employed is 

-0.0003 for wings without dihedral and 0.0005 for wings with dihedral. Also, the absolute 
value of the maximum incidence angle difference was determined to be less than 4° at 
98 percent semispan. At the next inboard station, 94 percent semispan, the absolute value 
of the difference was reduced to less than 1°. Hence, the error is highly localized and 
could be accounted for by extrapolation of information inboard of the tip in the layout of a 
model. Thus the sensitivity to Ng is essentially the same as before. Consequently, 
further calculations presented herein for wings with dihedral use values of Ng based on 
the initial sensitivity for wings without dihedral. 
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Three-dimensional comparisons .- Two comparisons with available mean-camber- 
surface solutions will be made. The comparisons are for a high-aspect-ratio sweptback 
and tapered wing with a uniform area loading at ^ = 1.0 and M^o = 0.90 and a lower 
aspect-ratio trapezoidal wing with a = 1.0, spanwise elliptic loading at ^ = 0.35, and 
Moo = 0.40. 

Figure 6 presents the predicted results from the present method for the sweptback 
wing and compares these results with those from references 1 and 10. A comparison of 
the three solutions indicates that they are all in generally good agreement with the excep- 
tion of the results at — ^ = 0.05. The surprising result is that the present method and 

b/2 

the modified Multhopp method (ref. 1) agree as well as they do at this span station because 
of the known differences that exist between them near the plane of symmetry. The reason 
for the larger disagreement between the present method and that of reference 10 near 

= 0 is not clear, but this disagreement may be caused by the different Nf. values 

b/2 

utilized by the two methods. Reference 10 effectively uses an infinite number since over 
each infinitesimal span strip across the wing the method locates a single quadrilateral 
vortex around the periphery of the enclosed area. This vortex extends from the leading 
edge to the trailing edge and includes segments of the edges as well. For a uniform area 
loading, the trailing leg parts of the quadrilateral vortices cancel with adjacent spanwise 
ones all across the wing. This leaves only the edge segments to contribute to the induced 
flow field. The present method utilizes a numerical rather than a graphical solution in 
order to provide a general capability; hence, N^. values are limited as discussed pre- 
viously. Also, vortices are not placed around the leading and trailing edges in the pres- 
ent method. 

A comparison of the present design method with that of reference 1 is shown in fig- 
ure 7. The wing and loadings are the same as those used in figure 4. The local slopes 
and elevations determined by the two methods are in reasonably close agreement at the 
three spanwise locations detailed; however, an oscillatory trend is evident in the local 
slopes obtained from the method of reference 1 (fig. 7(a)). These oscillations apparently 
originate in the truncated sine series used in reference 1 to represent a uniform chord- 
wise distribution. Integration of the local slopes to obtain local elevations tends to sup- 
press the oscillations (fig. 7(b)); however, the local pressures depend upon the slope 
rather than the elevation. Consequently, the measured chordwise pressure distribution 
will demonstrate the same oscillatory character. A model built according to the design 
of reference 1 was tested (ref. 11), and the measured pressure distributions for a typical 
spanwise location (fig. 7(c)) indicate that indeed the oscillations are present. Presumably, 
similar measurements on a model designed by the present method would not behave in this 
manner since the input loadings are truly linear. 
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Force tests (ref. 12) of an essentially identical model indicate that the measured 

drag polar was tangent to Cj) = that is, the vortex drag was indeed a mini- 
mum at the design (or 100 percent leading-edge suction was obtained). It is pre- 

sumed from the small differences in local slope between the present method and the 
method of reference 1 that a similar result would be obtained for a design by the present 
method. 


Application to a Wing-Canard Combination 

The present method has been demonstrated by optimizing a wing-canard combination 
(fig. 8). The effects of varying the vertical separation and the moment trim point on the 
resulting drag, span loading, and mean camber surfaces are also illustrated. All surfaces 
are designed for C^ ^ = 0.2, a^ = 0.6, a^ = 0.8, and = 0.30 and have Cm = 0 
about the moment trim point. Figure 8 shows that for all vertical separations, moving the 
moment trim point forward increases the vortex drag over some range, and furthermore, 
increasing the out-of-plane vertical separation reduces the vortex drag. Of course, not 
all moment trim points utilized will produce a stable configuration. These variations 
illustrate the importance of balancing the lift between the two lifting surfaces so that for 
some reasonable moment trim point and vertical separation, the vortex drag will be at 
a minimum. The minimum point on each vortex drag curve occurs with the pitching- 
moment constraint not affecting the extremization. 

The idea of lift balancing is an interesting one and is explored further for a moment 

\ Y 

trim point corresponding to = 0.1. Figure 9 shows the individual and total span load- 

b/2 

Zf. 

ings for the wing-canard configuration at = 0 for various values of a^ and a^. 

b/2 

From these figures there are three important observations to be made: (1) The individual 
span loadings change in the anticipated direction with the changing chord loadings in order 
to meet the same Cl and constraints; (2) the total span loading does not change; 

(3) consequently, the vortex drag of the configuration is constant, as would be anticipated 
from Munk's stagger theorem. 

Figure 10 presents the individual span loadings with increasing vertical separation 

I \ 

I — — < 0 above the wing plane with a« = 0.6 and a^ = 0.8. There are three observa- 
\b/2 / 

tions which can be made from these results for increasing vertical separation: (1) The 
individual span loadings tend to become more elliptical; (2) consequently, the vortex drag 
decreases; (3) the individual lift contributions show only a little sensitivity to separation 
distance once the canard is above the wing, when compared with the coplanar results. 
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Figure 11 shows the effect of moving the moment reference point on the spanwise 

distribution of wing and canard incidence angle for — = -0.169. The general result 

b/2 

shows that moving the moment reference point aft reduces the amount of incidence -angle 
nonuniformity required on each planform. This reduction is attributed to the change in 
loading on the canard required to meet the pitching-moment constraint. 

Figure 12 shows the effect of varying the vertical separation on the spanwise dis- 

A V 

tribution of wing and canard incidence angles for — = 0.1. As expected, with increas- 

b/2 

1 ing vertical separation the incidence requirements on each planform are generally reduced 

Zp Zp 

and should tend to the free-air result as — ^ - -». Note that for = 0, the wing is 

b/2 b/2 

required to have severe incidence gradients near the canard tip at — ^ = 0.673. This 

b/2 

unrealistic result occurs because the canard tip vortex intersects the wing, thereby induc- 
I ing a strong downwash field inboard and a strong upwash field outboard. These large inci- 
I dence gradients indicate that large out-of-plane displacements are called for in this solu- 

1 tion. The preceding results are, however, academic and occur as a result of the planar 

wake assumption and do not account for any real-wing effects or canard-wake rollup. 

Zc 

Two additional canard positions were examined: one at = -0.0845 and the 

b/2 

Zf. 

other at — ^ = 0 with 20® of dihedral. In each position, as could be expected, the large 
b/2 

changes in incidence on the wing, which occur near the canard tip spanwise location, are 

Zc 

significantly reduced and approach those of the other solutions. This helps to con- 

b/2 

Zp 

firm that the earlier solution for = 0 is special, and the large incidence gradients 

b/2 

noted can be avoided by providing the canard with a small effective displacement relative 

A V Zp 

to the wing. Additional details of the solution with = 0.1 and = 0.0 are given 

b/2 b/2 

in appendix D in sample case 2. 

Figure 13 presents selected local elevations for the wing and canard designed in the 

A V Zp 

presence of one another and alone at = 0.1 and = -0.676. For the wing the pri- 

b/2 b/2 

mary effect of adding the canard is to increase the incidence angle of the wing to compen- 
sate for the canard downwash field. For the canard there is only a small effect of being 
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designed in the presence of the wing - a reduction (or increase) in the incidence required 
when the wing induced field is upwash (or downwash). When the surfaces were designed 
alone, the same individual Cl as obtained in the combination design was used, and the 

chord load fraction ^a^. or a^^ was retained. Thus, the only loading variable between 
the two sets was the span loading, which was kept elliptical for the planform alone designs. 

Application to Tandem Wing Design 

This design method has been employed in the determination of the local elevation 
surfaces for a tandem wing. Figure 14 shows a sketch of a tandem wing configuration 
and selected results taken from the wind-tunnel tests made with a model based on this 
design at a Mach number of 0.30 (ref. 13). At Cl = 0.35 the vortex drag increment 
is correctly estimated. The measured Cm is slightly positive (0.02). Reference 13 
states that a part of the Cm error ^Cm should be zeroj is a result of a difference in 
the fuselage length between the designed and constructed model. 

Design of a Wing-Winglet Configuration 

Figure 15 presents a wing-winglet combination of interest along with pertinent 
aerodynamic characteristics and local elevations obtained from the present method. For 
comparison these same items are calculated with a program modification that adds a root- 
bending-moment constraint to produce the same moment that would be obtained on the orig- 
inal wing extending to the plane of symmetry but without its basic wingtip. The assumed 
span loading is elliptical. (See appendix F for a discussion of the root-bending-moment 
constraint.) The force and moment coefficients are based on the wing outside of a repre- 
sentative fuselage and without the basic wingtip. 

The results of this comparison are as follows: (1) The root-bending-moment con- 
straint increases the vortex drag slightly because of the changes in the c^c distribution 
required; (2) the differences in local elevations are confined primarily to the outer 50 per- 
cent semispan and are mainly due to the differences in the incidence angles; (3) significant 
amounts of incidence are reqviired in the winglet region with or without the root-bending- 
moment constraint. 

Additional details of the solution without the root-bending-moment constraint are 
provided in appendix D in sample case 1. 

The local elevation surfaces for a wing having both an upper and lower winglet can 
also be designed with this program when the two-planform option is employed. However, 
for such a configuration it is recommended that the pitching-moment constraint be 
dropped. 
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CONCLUDING REMARKS 


A new subsonic method has been developed by which the mean camber (local eleva- 
tion) surface can be determined for trimmed noncoplanar planforms with minimum vortex 
drag. This method employs a vortex lattice and overcomes previous difficulties with 
chord loading specification. This method designs configurations to have their local mid- 
surface elevations determined to yield the span load for minimum vortex drag while 
simultaneously controlling the pitching- moment or root-bending-moment constraint at 
the design lift coefficient. This method can be used for planforms which (1) are isolated, 
(2) are in pairs, (3) include a winglet, or (4) employ variable sweep, but only at a speci- 
fied sweep position. 

Results obtained with this method are comparable with those from other methods 
for appropriate planforms. The versatility of the present method has been demonstrated 
by application to (1) isolated wings, (2) wing-canard configurations, (3) a tandem wing, and 
(4) a wing-winglet configuration. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
April 7, 1976 
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APPENDIX A 


VORTEX-LATTICE COMPUTER PROGRAM FOR DETERMINATION OF 
MEAN CAMBER SURFACE (LANGLEY COMPUTER PROGRAM A4062) 

Basic Concepts and Limitations 

The vortex-lattice method is used in this computer program to determine the 
mean camber surfaces of planforms at subsonic speeds. This method assumes steady, 
irrotational, inviscid, incompressible, attached flow. The effects of compressibility are 
represented by application of the Prandtl-Glauert similarity rule to modify the planform 
geometry. Potential flow theory in the form of the Biot-Savart law is used to represent 
disturbances created in the flow field by the lift distribution of the planform. Those ver- 
tical displacements which occur in the configuration as a result of either dihedral or non- 
coplanar planforms are taken into account in the implementation of the Biot-Savart law. 
However, local displacements above or below the chord line at any spanwise position are 
ignored in the implementation. 

The planform is divided into many elemental panels. Each panel is replaced by a 
horseshoe vortex. This horseshoe vortex has a vortex filament across the quarter-chord 
of the panel and two filaments streamwise, one on each side of the panel starting at the 
quarter-chord and trailing downstream in the free-stream direction to infinity. Figure 1 
shows a typical horseshoe -vortex representation of a planform. 

The lifting-surface planform is represented for the computer program by a series 
of up to 24 straight segments which are positioned counterclockwise around the perimeter 
of the left half of the planform. Lateral symmetry is presumed. The lines start on the 
leading edge at the most inboard y- location, go along the leading edge to the left tip of the 
planform, return along the trailing edge, and end on the trailing edge of the most inboard 
y- location. The preciseness of the x and y Cartesian coordinates and dihedral angles, 
given as input data, determines the accuracy of the planform representation. It is recom- 
mended that the planform coordinates listed in the second group of geometry output data 
given in appendix C be plotted and examined after each computation to verify the accuracy 
of the planform representation. This check should be made before using the aerodynamic 
or local elevation output data. 

There are a number of restrictions and limitations in the application of this com- 
puter program. These limitations are discussed in detail in the program description and 
are noted with the appropriate input variables in appendix B. For the convenience of the 
program user, a complete list of restrictions and limitations is presented here. 
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The restrictions in the first group apply to all planforms and are as follows: 

(1) A maximum of two planforms may be specified. For examples, see sample 
case 1 for one planform and sample case 2 for two planforms. 

(2) A maximum of 24 straight-line segments may be used to define the left half of a 
planform. The lateral separation of the ends of these lines can be critical when the horse- 
shoe vortices are laid out by the computer program. For details of the manner in which 
the program handles the lateral separation, see Part I, Sections 2 and 3 under "Program 
Description." 

(3) The maximum number of horseshoe vortices on the left side of the configuration 
plane of symmetry is 400. When two planforms are specified, the sum total of the vor- 
tices in both is limited to 400. Within this limit, the number of horsehsoe vortices in any 
chordwise row may vary from 1 to 20 and the total number of chordwise rows may vary 
from 1 to 50. For examples, see the sample cases in appendix D. 

The limitations that apply only to variable- sweep planforms are as follows: 

(1) There should always be a fixed- sweep panel between the root chord and the outboard 
variable -sweep panel; (2) the pivot cannot be canted from the vertical; (3) no provisions 
have been made for handling dihedral in the geometry calculations for the variable-sweep 
panel or at the intersection of this panel with the fixed position of the wing. Restrictions 
on allowed values or codes for individual items of input data are described in appendix B. 

The calculations presented herein were made with a computer which used approxi- 
mately 15 decimal digits. For other computers with fewer significant digits, it may be 
necessary to use double precision for some of the calculations. In addition, it may be 
necessary to change some of the tolerances used in the program. These tolerances are 
given in the program listing. 

Program Description 

This FORTRAN program is used to compute the local elevation shapes of multiple 
lifting planforms and is divided into three parts. Part I contains the geometric calcula- 
tions, Part n contains the circulation term calculations, and Part in contains the final 
output terms and answer listings. These three parts describe the three types of compu- 
tations performed in the FORTRAN computer program. The input data are described in 
detail in appendix B, and the output data are described in detail in appendix C. Two sam- 
ple cases are given to illustrate the use of the program. Listings of the input data and 
computed results for these sample cases are given in appendix D, and the FORTRAN com- 
puter program is given in appendix E. 
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Part I - Geometry Computation 

The first part of the program is used to compute the geometric arrangement 
required to represent the planform by a system of horseshoe vortices and is divided into 
three sections. In Section 1, a description of the planform (group one of the input data in 
appendix B) is read into the computer. In Section 2, configuration details (group two of 
the input data) are read into the computer. In Section 3, the horseshoe vortex lattice is 
laid out. When two planforms are used to describe a wing-body-tail configuration, each 
of these sections is repeated for the second planform. At the beginning of the geometry 
computation, a data card is read which describes the number of planforms (either 1 or 2), 
the number of configurations for which values are to be computed, and the reference val- 
ues for chord and area. 

Section 1 - Reference Planform: 

The planform is described by a series of straight lines which are projected onto the 
X-Y plane from the deflected planform as shown in figure 1 for a double-delta type plan- 
form. The primary geometric data are the locations of the intersections of the perim- 
eter lines, the dihedral angles, and an indication as to whether the lines are on a fixed or 
movable panel. (See ref. 4 for an example.) The pivot location is also required for a 
variable-sweep planform. These data are described in group one of the input data (appen- 
dix B). For variable-sweep wings, the planform used for input should be the configuration 
with the movable panel in a position where the maximum number of lines required to form 
its perimeter is exposed. 

Section 2 - Configuration Computations: 

The particular configuration for which the local elevation surface is sought is 
described by group two input data which are read in this section. These data include the 
following quantities: an appropriate configuration number, the number of horseshoe vor- 
tices chordwise, the nominal number of chordwise rows of vortices spanwise, the Mach 
number, the particular lift coefficient at which the local elevation surface is desired, and 
the sweep angle of the outboard panel for variable -sweep wings. 

The number of horseshoe vortices used in each chordwise row (SCW) must be con- 
stant across the span. Simply indicate the number on the configuration card and this 
value will be used on each planform of the group one input. For all but the most simple 
planforms, the program adds some extra rows of horseshoe vortices. (This is discussed 
in Section 3.) As a result, the number of chordwise rows actually laid out (SSW) is usually 
greater than the nominal number of rows (VIC), and it takes one complete run through the 
program to determine the exact number and location of the rows. If variations in the 
basic wing planform are desired for additional computer cases, the entire computer pro- 
gram must be rerun with all geometry data and the appropriate changes in any of the 
aforementioned variables in the group two input data. 
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For a variable -sweep planform, the angle which describes the sweep should be on 
the leading edge of the movable panel adjacent to the fixed portion. The intersection 
points and sweep for the planform in the desired position are then computed. For a fixed 
planform, the sweep-angle specification is not required because the program will use the 
vmaltered basic planform. The planform breakpoints are checked to see whether the 

spacing between any consecutive pair in the spanwise direction is less than If this 

occurs, the points are adjusted to coincide with each other. The adjustment is necessary 
to avoid a poorly conditioned matrix which could result in biased results for the 9z/ 8x 
terms. Although this adjustment is usually adequate for planforms with no dihedral, it 
may not be sufficient for a particular configuration with dihedral or for use of this pro- 
gram in computers which have fewer than 15 significant decimal digits. This problem is 
discussed in detail in Section 3. 

When two planforms are specified, the program compares the spanwise location of 
the breakpoints on both planforms inboard of the tip of the planform with the shorter semi- 
span. If all the breakpoints coincide spanwise, no action is taken. However, if one plan- 
form has a breakpoint which does not occur on the other planform, an additional brealq)oint 
is added to the other planform on its leading edge. This is done to force all trailing legs 
from the horseshoe vortices to occur at the same spanwise location, which keeps a trail- 
ing leg from one planform from passing too close to a slope point on the other planform 
and prevents unrealistic induced velocities at that slope point. 

The program determines the planform area and span projected to the X-Y plane and 
uses these values to compute the average chord. Planforms which have a constant angle 
of dihedral from the root chord to the tip chord have an average chord which is indepen- 
dent of dihedral angle. However, wings with more than one dihedral angle have an aver- 
age chord which is dependent on the individual dihedral angles. 

Section 3 — Horseshoe Vortex Lattice: 

In this section, the procedure by which the horseshoe vortex lattice is laid out is 
described. The planform is divided chordwise and spanwise along the surface into trap- 
ezoidally shaped elemental panels; one horseshoe vortex is assigned to represent each 
panel. The horseshoe vortices are the same as those described in reference 4 and one 
is sketched in figure 2 for a typical panel. The horseshoe vortex is composed of three 
vortex lines: a bound vortex which is swept to coincide with the elemental-panel quarter- 
chord sweep angle in the plane of the wing and two trailing vortices which extend chord- 
wise parallel to the free stream to infinity behind the wing. Figure 1 shows a typical 
chordwise row of horseshoe vortices on an arbitrary planform. The nominal width of 
these horseshoe vortices is the total semispan in the plane of the wing divided by the 
variable VIC. (See appendix B.) 
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The procedure for laying out the elemental panels and, consequently, the horseshoe 
vortices is to begin at the left leading-edge tip with a chordwise row of horseshoe vortices 
and then to proceed inboard toward the most inboard y- location of the wing. The actual 
spanwise locations of the chordwise rows of horseshoe vortices are adjusted so that there 
is always a trailing vortex filament at points where there are intersections of perimeter 
lines or breakpoints on the planform. This adjustment may cause the horseshoe vortex 
width to be narrower or wider than the nominal width. When a horseshoe vortex has one 
trailing vortex filament which coincides with a breakpoint, the width of the horseshoe vor- 
tex may vary from 0.5 to 1.5 times the nominal width. When both trailing legs coincide 
with breakpoints, the width may vary from a maximum of 1.5 times the nominal width to 

/rt 

a minimum width of — — , as described previously in Section 2, The number of chord- 

2000 

wise rows actually laid out is given by the variable SSW. 

In the chordwise direction, the horseshoe vortices are distributed uniformly and the 
number of vortices is given by the variable SCW. The maximum number of horseshoe 
vortices in the chordwise direction is 20, and in the spanwise direction the maximum total 
number of chordwise rows is 50 on a semispan. However, the total number of horseshoe 
vortices (the product of SCW and SSW) permitted by the program is 400 on the left half of 
a configuration. The exact number generated by the program depends on the value of VIC 
and SCW and on the details of the planform. As many as one additional chordwise row of 
horseshoe vortices may be generated by the program at each breakpoint outboard of the 
root. Wings with dihedral must always have at least two horseshoe vortices chordwise; 
wings without dihedral may have only one. 

The Prandtl-Glauert correction factor is applied to the x- coordinates and the tan- 
gents of the sweep angle of the horseshoe vortices at this point to account for compress- 
ibility effects. 

Part II - Vortex Strength Computation 

The vortex lattice laid out in Part I is not employed to determine the vortex 
strengths, but instead is utilized to find the local elevation shapes (Part III) because of 
the smaller computer resource requirements. (See the section "Lift, Pitching-Moment, 
and Drag Contributions" for additional discussion.) The solution for the vortex strengths 
is accomplished in the Trefftz plane by using the one or two lifting lines which may be 
bent. These lines are 'divided into equal segments, with 50 divisions per planform semi- 
span used for the planform with the larger true length. In case of two planforms of 
unequal length, the number of equal segments assigned to the shorter lifting line is pro- 
portional to the length ratio of the two planforms. These segments are laid out from in- 
board to outboard on the lifting lines. For the shorter lifting line, a small portion near 
the tip may not be included but will always be less than 2 percent of the larger semispan 
true length because of the use of whole equal segments. 
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After the optimization is performed, in which the spanwise scaling factors are 
determined based on a Trefftz plane solution, these scaling factors are interpolated back 
to the original spanwise paneling layout for the vortex lattice. It is these results which 
provide the multipliers for the chordwise shapes. They lead to the computation of the 
span loadings. Cl, and Cm developed for each planform. The circulations are listed 
and then employed in Part III. 

Part III - Local Elevation Shape Computation 

The vortex strengths determined in Part II and the influence coefficient matrix 
based on the original paneling (see fig. 1) are used in this part of the program to com- 
pute the local slope at the midspan three-quarter-chord location of each elemental panel 
(called the slope point in fig. 2) by employing equation (1). By using cubic splines to 
interpolate between the local slopes, the local elevation shape at each spanwise location 
is determined by equation (28). Outside the range of slope points a constant extrapolation 
procedure is used to determine the integrand of equation (28). (See the section "Preci- 
sion” for a discussion of the extrapolation methods examined.) 
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INPUT DATA 
Group One 

The input data required for the reference planform are described in the order that 
they are called for by the computer program. All coordinates and sweeps should be given 
for the left half of the wing planform. The axis system used is given in figure 1 and any 
consistent set of units is acceptable. (The output will be in terms of the input units.) 

The X-axis coincides with the plane of symmetry and is positive pointing into the wind; 
the Y-axis is positive pointing along the right wing. The origin of the axis system may 
lie anywhere along the plane of symmetry and determines the trim point for the two- 
planform solution. All the cards use a format of 8F10.6 for group one data. 

Data on the first card are for the four named variables and are to be supplied in the 
following order: 

PLAN number of planforms for the configuration; use 1 or 2; this sets the 

maximum number for the IT variable used subsequently 

TOTAL use 1 for this field 

CREF reference chord of the configuration; this chord is used only to nondi- 

mensionalize the pitching-moment terms and must be greater than 
zero 

SREF reference area of the configuration; this area is used only to nondimen- 

sionalize the lift, drag, pitching moment, and root bending moment 
and must be greater than zero 

The data required to define each planform are then provided by a set of cards. The 
initial card in this set is composed of the following data: 

AAN(IT) number of line segments used to define left half of a wing planform 

(does not include root chord); a maximum of 24 line segments may 
be used 

XS(IT) x-location of the pivot; use 0 on a fixed wing; the axis system used is 

given in figure 1 
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YS(IT) y-location of the pivot; use 0 on a fixed wing 

RTCDHT(IT) vertical distance of particular planform being read in with respect to 
the wing-root-chord height; use 0 for a wing 


The rest of this set of data requires one card for each line segment used to define 
the basic planform (variable AAN(IT)). All data described below are required on all 
except the last card of this set; the last card uses only the first two variables in the fol- 
lowing list: 

XREG(I,IT) x-location of ith breakpoint; the first breakpoint is located at the inter- 

section of the left wing leading edge with the root chord; the break- 
points are numbered in increasing order for each intersection of 
lines in a counterclockwise direction 


YREG(I,IT) y-location of ith breakpoint 

DIH(I,IT) dihedral angle (degrees) in Y-Z plane of line from breakpoint i to 

i + 1, positive upward; along a streamwise line, the dihedral angle is 
not defined; use 0 for these lines; the dihedral angle will have the 
same sign and magnitude along the leading and trailing edges of a 
planform over the same spanwise extent 


AMCD move code; this number indicates whether the line segment is on the 

movable panel of a variable- sweep wing; use 1 for a line which is 
fixed or 2 for a line which is movable 


Group Two 

Two sections of data form the group two data. The first section is a single card 
which describes the details of the particular configuration for which the mean camber 
surface is desired. This card requires a format of 5F5.1, 2F10.4. The second section 
is used to supply the fractional chordwise locations where the chord load changes from 
a constant value to a linearly varying value toward zero. This card uses a format of 
8F10.4. 

Section one data are to be supplied in the following order: 

CONFIG arbitrary configuration number which may include up to four digits 
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sew number of chordwise horseshoe vortices to be used to represent the 

wing; a maximum value of 20 may be used; do not set to zero 

VIC nominal number of spanwise rows at which chordwise horseshoe vorti- 

ces will be located; the variable VIC must not cause more than 50 
chordwise rows of vortices to be used by the program to describe 
the left half of the configuration; in addition, the product of SSW and 
sew cannot exceed 400; the use of the variable VIC is discussed in 
detail in Part I, Section 3 of appendix A 

MACH Mach number; use a value other than 0 only if the Prandtl-Glauert com- 

^ \ 

pressibility correction factor ^/3 = ||l - Moo^j is to be applied; it 
should be less than the critical Mach number 

CLDES design lift coefficient for lifting system 

SA(1) variable sweep angle of the first planform; specify leading-edge sweep 

angle (degrees) for the first movable line adjacent to the fixed portion 
of the planform; for a fixed planform this quantity may be omitted 

SA(2) variable sweep angle for the second planform 

Section two data consist of two quantities: 

XCFW fractional chord location where the chord load changes from a constant 

value to a linearly varying value toward zero at the trailing edge of 
the first planform; this is the same as the symbol a used in the 
body of the paper 

XCFT fractional chord location where the chord load changes from a constant 

value to a linearly varying value toward zero at the trailing edge of 
the second planform; this is the same as the symbol a used in the 
body of the paper; if only one planform is present, the variable XCFT 
should be omitted from the input data 

Guidelines for Program Use 

The following guidelines for the use of this program have been developed from iso- 
lated wing studies using the solution technique for configurations without dihedral: 
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(1) More than 10 and perhaps as many as 20 horseshoe vortices are needed along a 
chord to assure a good solution for the mean camber surface. 

(2) At least 10 chordwise rows of horseshoe vortices should be used along a semi- 
span. More chordwise rows can be used to save interpolating time, although they will not 
necessarily yield a better solution. 

(3) Uniform spacing of horseshoe vortices chordwise and of the chordwise rows 
spanwise is preferred. 
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OUTPUT DATA 

The printed results of this computer program appear in three sections: geometry 
data, aerodynamic data, and local elevation data. 

Geometric and Aerodynamic Data 

The geometry data are described in the order that they are found on the printout. 
The first group of data describes the basic planform, stating the numbers of lines used 
to describe the planform, the root-chord height, and the pivot position and then listing 
the breakpoints, sweep and dihedral angles, and move codes. These data are a listing 
of the input data except for the sweep angle, which is computed from the input data. 

The second group of data describes the particular planform for which the local 
elevation data are being computed. Included are the configuration number, the sweep 
position, a listing of the breakpoints of the wing planform in terms of (x,y,z), the sweep 
and dihedral angles, and the move codes. These data are listed primarily for variable- 
sweep wings to provide a definition of the planform where the outer panel sweep is dif- 
ferent from that of the reference planform. 

The spanwise scale factors and the term 

Normal induced velocity 

(Free-stream velocity) [cos (local dihedral angle)] 

are listed between the second and third groups of data if the configuration has dihedral. 

The third group of data presents a detailed description of the horseshoe vortices 
used to represent the planform. These data are listed in eight columns, with each line 
describing one elemental panel of the wing. The following items of data are presented 
for each elemental panel: 

X C/4 x-location of quarter-chord at horseshoe vortex midspan 

X 3C/4 x-location of three-quarter-chord at horseshoe vortex 

midspan; this is the x-location of the slope point 

Y y-location of horseshoe vortex midspan 

Z z -location of horseshoe vortex midspan 
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S semiwidth of horseshoe vortex 

C/4 SWEEP ANGLE sweep angle of quarter-chord 

DIHEDRAL ANGLE dihedral angle of elemental panel 

GAMMA/U AT CLDES = P/U distribution at the design 

The fourth group of data presents the following geometric data: 

REF. CHORD reference chord of wing 

C AVERAGE average chord (true planform area divided by true span) 

TRUE AREA true area computed from planform listed in second group 

of geometry data 

REF. AREA reference area 

B/2 largest true semispan of the planforms listed in second 

group of geometry data 

REF. AR reference aspect ratio computed from reference planform 

area and true span 

TRUE AR true aspect ratio computed from true planform area and 

true span 

MACH NUMBER Mach number 

The following aerodynamic data are given: 

CL*C c^c, span loading 

CL DESIGN Cl design Cl 

CL COMPUTED total Cl actually developed from the interpolated span- 

wise scaling results 


39 



APPENDIX C 


CM COMPUTED total Cj^ actually developed from the interpolated span- 

wise scaling results 

CD V Cjj vortex drag coefficient based on the far-field solu- 

tion at ^ 

Local Elevation Data 

This section contains the local elevation solutions along the semispan of up to two 
planforms. An explanation of the variables listed is as follows: 


Y 

y, physical spanwise location 

Y/B/2 

y 

, fraction physical spanwise location based on semi- 

b/2 

span of larger planform 

CHORD 

physical chord at y 

DZ/DX 

dz/dx, slope of local elevation curves along the chord 

x/c 

fractional chordwise distance measured from the leading 
edge, positive aft 

z/c 

z/c, local elevation normalized by the chord measured 
with respect to the local trailing edge, positive down 

DELTA X 

(x/c) (chord) 

DELTA Z 

(z/c) (chord) 
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SAMPLE CASES 

Input data, sketches, and output data for two sample cases are presented in the fol- 
lowing order: 


Sample case 

Configuration 

Item 

Page 

1 

100 

Input data 

42 



Sketch 

42 



Output data 

43 

2 

2 

Input data 

70 



Sketch 

70 



Output data 

71 


Using the same solution technique leads to the central processing time for a config- 
uration generally increasing as the square of the increase in the number of horseshoe 
vortices used to represent the left half of the planform. Some typical times for the 
sample cases with a Control Data Corporation 6600 computer system are as follows: 


Sample case 

Solution technique for 
configuration having — 

Number of 
horseshoe vortices 

CPU time, sec 

1 

Dihedral 

340 

140 

2 

No dihedral 

400 

183 
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Input Data and Sketch for Sample Case 1 


COLUMN NUMBERS FOR INPUT DATA 
0000000001 11 1 1 1 1 1 1 1 ???2222???33 13131 ‘^3344AA4AAA44S5SS5SSS556f»666^6hA677 7777 7777A 
1 234567990 12145678Q0 1234567890 1 ?345678Rn 1234567PQ01 2 145^S7A9 o1234S67«901?34S^7R90 

GROUP ONE DATA 


1. 

1. 

18.145 

1782.272 

8. 

0. 

0. 

0. 

26.68 

-0. 

6 . 

1. 

-20.52 

-60. 

77.5 

1. 

-22.82 

-60.65 

77.5 

1. 

-29.06 

-61.861 

0. 

1 . 

-30.58 

-61.861 

77.5 

1 . 

-27.72 

-60.65 

77.5 

1 . 

-27.54 

-60.0 

6 . 

1 . 

-12.12 

-24. 

6 . 

1 * 

-7.92 

-0. 




GROUP TWO 

DATA 


100. 

20. 18. 0 

.8 0.5 



1.0 


No root -bending-moment constraint is employed. 
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APPENDIX D 


Input Data and Sketch for Sample Case 2 


COLUMN NUMBERS FOR INPUT DATA 
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c/4 SWEEP DIHEDRAL GAMMA/U AT 
ANGLE ANGLE CLOES* .2000 
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FORTRAN PROGRAM LISTING 

This program was written in FORTRAN IV language, version 2.3 for the Control 
Data Corporation series 6000 computer system with SCOPE 3.0 operating system and 
library tape. Minor modifications may be required prior to use on other computers. 

The program is written using UPDATE and PROGRAM stepping. These features allow 
the program storage requirements to vary from 51000g to 1120008 words, depending on 
the matrix conditioning and the solution technique for the aerodynamic characteristics. 
The solution technique for configurations without dihedral uses PROGRAM CIRCULl and 
51000g words; the solution technique for configurations with dihedral uses PROGRAM 
CIRCUL2 and 63000g words for a well-conditioned matrix and uses PROGRAM CIRCUL3 
and 1120008 words for an ill-conditioned matrix. The selection takes place automatically 
and is dependent on the geometry of the configuration and the vortex-lattice layout. 

This computer program consists of four basic PROGRAM steps, three OVERLAYS 
and seven SUBROUTINES. Each PROGRAM, OVERLAY, and SUBROUTINE is identified 
in columns 73 to 75 by a three-letter abbreviation. In addition, each of these parts is 
sequenced with a three-digit number in columns 77 to 79. The following table is an index 
to the program listing: 


Name of part 

Abbreviation 

Page 

PROGRAM GEOMTRY 

GEO 

109 

OVERLAY 0 (WINGTL) 
PROGRAM WINGAL 

1 DGO 

119 

SUBROUTINE FTLUP 

TLU 

120 

SUBROUTINE SIMEQ 

SEQ 

122 

SUBROUTINE DRAGSUB 

DGS 

124 

OVERLAY 1 (WINGTL) 

DGl 

125 

PROGRAM CIRCULl 



OVERLAY 1 (WINGTL) 

\ DG2 

130 

PROGRAM CIRCUL2 
OVERLAY 1 (WINGTL) 

J 

1 DG3 

135 

PROGRAM CIRCUL3 
SUBROUTINE GIASOS 

) 

GIA 

140 

OVERLAY 2 (WINGTL) 

> zoc 

148 

PROGRAM ZOCDETM 
SUBROUTINE INFSUB 

J 

INF 

150 

SUBROUTINE SPLINE 

SPL 

151 

SUBROUTINE TRIMAT 

TRI 

153 

PROGRAM DUMMY^ 

DUM 

153 


^The PROGRAM DUMMY is for default purposes of 
PROGRAM GEOMTRY. 
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J0tj,l» 1000*063000«1000. R4310 lOOllO 

JSER. LAMAR* JOHN £ O00503AOON 38S10 

'<ORFL. 

jPDATE(F.I*N,CtL=0) 

REi<IND(NEWPL) 

JPDATE (Q*P=NEWPL*C.L=0) 

-fU'^(S*»«COMPILt) 

SETINDF. 

lGO. 

RE*(IMO(NEWPL) 

REi*IND (TAPE50) 

JHDAfE (Q, I=T APE50*P=NEWPL.L=0) 

RUN(S.**COMPILt.*6LO) 

SETINDF. 

3LO. 

EXIT. 
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•DECK 


VLHCGEOM 

program GEOMTRY ( 1 NPUT » OUTPUT. T APE5= I NPjTtTAPE6=OUTPUT»TAPE?5»TAPE5GEO 1 


10) 


geo 

2 

OIMENSIOM X«EF(25). YHEF(25). SAR(25). A(25). »SAR(25). X(2S). 

Y(2GE0 

3 

15), B0TSV(2), SA(2). VB0R0(51), SPY(50,2), KFX(2). 

lYL (50,2) , 

IYT(GE0 


250,2) 


GEO 

5 

COMMON /ALL/ SOT ,M , BET A , PTEST ,0TEST , TBlSCW ( 50 ) , 0 ( *,00 ) , PN (-.OO ) , 

PV(^GE0 

6 

100) ,S(A00) ,PSI (400) ,PH1 (50) ,ZH(50) ,NSSi^ 


GEO 

7 

COMMON /MAI NONE/ I COOEOF, TOT AL , A AN (2 ) ,XS(2) , YS (2) ,KFCTS(2) ,XRE6(25GE0 

S 

1,2) ,YREG(25,2) , A«EG(25,2) ,0IH(25,2) ,MC0(25,2) ,XX (25 

,2) ,YY (25,2) ,ASGE0 

9 

2(25,2) ,T Two (25,2) ,MMC0(25,2) , AN ( 2) ,ZZ ( 25, 2 ) , IFL AG 


GEO 

10 

COMMON /ONEThRE/ twist (2) ,cref,ssef,cave,cldes,strue,ah, arthue 

.rtcgeo 

11 

lOHT (2) ,C0 nFIG,NSSWSV(2) ,MSV(2) ,KB0T,PLAN, IPLAN.MACH 

,SSWWA(50) 

geo 

12 

COMMON /CCRRDO/ CHORD ( 50 ) , XTE ( 50 ) ,KB I T ,TSP AN , TSP AN A 


GEO 

13 

REAL MACH 


geo 

14 

REWIND 50 


GEO 

15 



GEO 

17 

PART ONE - geometry COMPUTATION 


GEO 

U 



GEO 

19 

SECTION ONE - INPUT OF REFERENCE WING 

position 

GEO 

20 



geo 

21 



GEO 

22 

icooeof=o 


GEO 

23 

TOTAL=PT£ST=QTEST=TWIST (1)=TWIST(2)=0. 


GEO 

24 

IF (TOTAL. EQ. 0. ) R TCDHT ( 1 ) =RTCDHT ( 2 ) =0 .0 


GEO 

25 

YTOL=l .E-10 


GEO 

26 

AZY=1.E*13 


GEO 

27 

PIT=1 .57079b3 


GEO 

28 

RAO=57. 29578 


geo 

29 

IF (TOTAL. GT.O.) GO TO 7 


GEO 

30 



GEO 

31 



geo 

32 

SET PLAN EQUAL TO 1. FOR A WING ALONE COMPUTAION - 

EVEN FOR A 

GEO 

13 

VARIABLE SWEEP WING 


geo 

34 

SET plan equal to 2. FOR A WING - TAIL COMBINATION 


geo 

35 



GEO 

36 

SET TOTAL EQUAL TO ThE NUMBER OF SETS 


GEO 

37 

OF GROUP TWO DATA PROVIDED 


GEO 

18 



GEO 

39 

READ (5,98) PLAN, TOT AL.CREF.SREF 


GEO 

40 

IF(ENOFIlE 5) 93,1 


geo 

41 

1PLAN=PLAN 


geo 

42 



geo 

43 



geo 

44 

SET AAN(IT) equal to the maximum number of curves required to 

geo 

^♦5 

DEFINE THE PlAnFORM PERIMETER OF THE (IT) plANFORM. 


GEO 

46 



GEO 

47 

SET RTCOHT(IT) equal TO THE ROOT CHORD HEIGHT OF THE LIFTING 

geo 

48 

SURFACE (IT), whose PERIMETER POINTS ARE BEING HEAD 

IN, with 

GEO 

49 

RESPECT TO The WING ROOT CHORD HEIGHT 


GEO 

50 



geo 

51 

WRITE (6,96) 


geo 

52 

DO 6 IT=1,IPLAN 


geo 

53 

READ (5,98) A AN ( I T ) , XS ( I T ) , YS ( I T ) ,RTCDHT ( I T ) 


GEO 

54 

N=AAN (IT) 


GEO 

55 

Nl=N* 1 


geo 

56 

MAK = 0 


GEO 

57 

IF (IPLAN.EQ.l) PRTCON=10H 


GEO 

58 

IF (IPLAN.EQ.2.AND.IT.E0.1) °RTCON=10H FIRST 


GEO 

59 

IF (IPLAN.E0.2.AND.IT.E0.2) PRTC0N=10H SECOND 


GEO 

60 
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2 


3 


9 


10 


11 


12 


»(RITE (6»97) PRTCON»NtRTCDHT ( IT) *XS ( IT) ♦ YS( IT) 


GEO 

61 

WRITE (6*109) 


GEO 

62 

DO 5 I=1*N1 


GEO 

63 

READ (5*98) XREG( I ♦ IT) *YREG( I ♦ IT ) «DIH( I . IT) • AMCO 


GEO 

64 

MCD(I*IT)=AMCD 


GEO 

65 

IF (I.EO.l) GO TO 5 


GEO 

66 

IF (^(AK.NE.O.OR.MCD(I-1*IT) .NiE.a) go to i 


GEO 

67 

MA<=I-1 


GEO 

68 

IF (AB5(YREG(I-1*1T)-YREG(I*IT) ) .LT.YTDL) GO TO 3 


GEO 

69 

AREG(I-1*IT) = (XREG(1-1.IT)-XREG( NIT) ) / ( YREG ( l-l * I T ) -YREG ( I 

f IT) ) 

GEO 

70 

ASWP = ATAN ( AREG( I-l * IT) ) *RAD 


GEO 

71 

GO TO A 


GEO 

72 

YRE6(I*IT)=YREG(I-1*IT) 


GEO 

73 

AREG(I-NIT)=A7Y 


geo 

74 

ASWP=90. 


GEO 

75 

J=I-1 


geo 

76 



GEO 

77 

write PUANFORM perimeter points and angles 


GEO 

78 



GEO 

79 

WRITE (6*106) I*XkE 6( J* lT)*YSEG(J*IT)*ASwp*DlH(J*in *MCU( J. 

IT) 

GEO 

HO 

DIH(J*IT)=TAN(niH( J*IT)/RAD) 


GEO 

81 

CONTINUE 


geo 

82 

KFCTS( IT)=MAK 


geo 

83 

WRITE (6*106) N1*XREG(N1.IT).YREG(N1*IT) 


geo 

84 

CONTINUE 


geo 

85 



GEO 

86 

PART 1 - SECTION i 


GEO 

87 

READ GROUP 2 DATA AND COMPjTE DESIRED wING POSITION 


geo 

88 



GEO 

89 



GEO 

90 

SET SA(1)*SA(2) EfJJAL TO ThE SWEEP ANGLE* IN DEGREES. FOR ThE FIRSTGEO 

91 

CURVE (S) THAT CAN CHANGE SWEEP FOR EACH PLANFORM 


geo 

92 



geo 

93 

READ ( 5. 1 OS) CONE IG*SCW*VIC*MACH* cedes *SA ( 1 ) ,SA (2) 


GEO 

94 



GEO 

95 

WRITE (6*99) CONFIG 


GEO 

96 

IF(ENDFIlE 5) 93.8 


GEO 

97 

IF (PTEST.NE.O..ANO.QTEST.NE.O.) GO TO 95 


GEO 

98 

IF (SCW.EO.O.) GO TO 10 


GEO 

99 

DO 9 1=1.50 


GEO 

100 

TBLSCW(D=SCW 


geo 

101 

GO TO 11 


geo 

102 

READ (5.98) STA 


geo 

103 

NSTA=STA 


geo 

104 

read (5.98) (TBLSCW( I ) .TBUSCW( 1*1 ) .TBLSCW ( 1*2) .TBLSCW( N3) . 

TBLSCW (GEO 

105 

1I*A) .TBLSCW ( 1*5) .T3LSCW( 1*6) .TBLSCW (1*7) . 1 = 1 .NSTA.8) 


GEO 

106 

DO 37 IT=1.IPLAN 


geo 

107 

N=AAN( IT) 


GEO 

108 

N1=N*1 


GEO 

109 

00 12 1=1.N 


geo 

110 

XREF(I)=XREG(I.IT) 


GEO 

111 

YREF(I)=YREG(I.IT) 


geo 

112 

A(I)=AREG(1.IT) 


GEO 

113 

RSAR(I)=ATAN(A(1) ) 


GEO 

114 

IF ( A( 1 ) .EQ.AZY) RSAR(I)=PIT 


GEO 

115 

CONTINUE 


GEO 

116 

XREF(N1)=XREG(N1*IT) 


GEO 

117 

YREF(N1)=YREG(N1.IT) 


GEO 

118 

IF (KFCTS( IT) .GT.O) GO TO 13 


GEO 

119 

X = 1 


GEO 

120 


no 
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13 

14 


C 

c 


15 

16 


17 


18 

C 

c 


c 

19 


20 

C 

c 

21 

22 


23 


24 


SA(IT)=WSA.<(l)»RAO 


geo 

121 

60 TO 14 


geo 

122 

K=KFCTS( IT) 


geo 

123 

WRITE (6*102) K*SA(IT)«IT 


geo 

124 

SB=SA( IT)/RAD 


GEO 

12S 

IF ( AflS (SH-RSAR <K) ) .6T. (, 1/RAO) ) 60 

TO 17 

geo 

126 

REFERENCE PLANFORM COORDINATES ARE ! 

STORED UNCHANGED FOR WINGS 

GEO 

127 

without CHANGE IN SWEEP 


GEO 

12B 

DO 16 1=1 *N 


GEO 

129 

X(I)=XREF(I) 


GEO 

130 

Y<1)=YREF(I) 


GEO 

131 

IF (RSAR(I) .EQ.PIT) 60 TO 15 


GEO 

132 

A(I)=TAN(RSAR(D) 


GEO 

133 

GO TO 16 


geo 

134 

A( I ) =AZY 


GEO 

135 

SAR( I ) =RSAR < I ) 


geo 

136 

X (Ml ) =X»EF (N1 ) 


GEO 

137 

Y (N1 ) =YR£F (N1 ) 


GEO 

138 

60 TO 35 


GEO 

139 



GEO 

140 

CHANGES IN WING SWEEP ARE MADE HERE 


geo 

141 



GEO 

142 

IF (MCD(X. IT) .NE.2) GO TO 94 


geo 

143 

KA=K-1 


GEO 

144 

<[ 

II 

00 

o 

Q 


GEO 

14S 

X (I ) =XREF ( I ) 


GEO 

146 

Y(I)=YREF(I) 


GEO 

147 

SAR ( I ) sRSAR ( I ) 


GEO 

148 

determine leading edge intersection 

BETWEEN FIXED AND VARIABLE 

GEO 

149 

Sweep wing sections 


geo 

ISO 

SAR (K ) =Sd 


GEO 

ISl 

A(K)=TAN(Sb) 


GEO 

152 

SAI=SB-RSAR(K) 


GEO 

153 

X(K*1 ) =XS (IT) ♦ (XREE (K*l )-XS( IT) ) *C0S(SAI ) ♦ ( YREF ( K* 1 ) -YS ( I T ) ) »SIN (SGEO 

1S4 

lAl) 


geo 

155 

y(k*1)=yS(it)* (yrep(k*1)-ys(it) ) »cos(sai )- (xref(k*d-xs( it > ) »sin(sgeo 

156 

lAI) 


GEO 

157 

IF (ABS(SH-SAR(K-1) ) .LT. (. 1/RAD) ) GO TO 19 

GEO 

158 

Y(K)=X(K*1)-X(K-1)-A(K)*Y(K»1)*A(K- 

1)*Y(K-1) 

GEO 

159 

Y (K ) =Y (K ) / ( A (K-1 )-A (K) ) 


GEO 

160 

X(K)=A(K)*X(K-1)-A(K-1)*X(K»1)*A(K- 

1)*A(K)*(Y(K»1)-Y(K-1) ) 

GEO 

161 

X(K)=X(K)/(A(K)-A(<-1) ) 


GEO 

162 

60 TO 20 


GEO 

163 

ELIMINATE EXTRANEOUS BREAKPOINTS 


geo 

164 

X(K)=XREF(K-1) 


GEO 

165 

Y(K)=YREF(K-1) 


GEO 

166 

SAR(K)=SAR(K-1) 


geo 

167 

K = K*1 


GEO 

168 

SWEEP The breakpoints on the variable sweep panel 

GEO 

169 

(IT ALSO KEEPS sweep ANGLES IN FIRST OR FOURTH QUADRANTS) 

GEO 

170 

K=K*1 


geo 

171 

SAR(K-1)=SAI*RSAR(<-1) 


GEO 

172 

IF (SAR(K-1 ) .LE.PIT) 60 TO 23 


geo 

173 

SAR (K-1 ) =SAR (K-1) -3, 1415927 


GEO 

174 

GO TO 22 


GEO 

175 

IF (SAR (K-1 ) ,0E. (-PIT) ) 60 TO 24 


geo 

176 

SAR ( K- 1 )=SAR( K-1) ♦3.1415927 


GEO 

177 

60 TO 23 


geo 

178 

IF ( (SAR(K-l) ),LT..O) GO TO 25 


GEO 

179 

IF (SAR(K-I)-PIT) 28*26.26 


GEO 

180 
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25 

26 

27 

28 

29 

30 


C 

C 

31 


33 

C 

C 


39 

35 


36 


37 


38 

39 
90 


IF (SAR(K-1)*PIT) 27,27.28 
A(K-1 )=AZY 
GO TO 29 
A (K-1 ) =-AZY 
GO TO 29 

A(K-1)=TAN(SAR(K-1) ) 

KK=MCD(K*IT) 

GO TO (31*30), KK 

Y (K)=YS ( IT) ♦ (YREF(<)-YS{ IT) )*COS (SAI ) - ( XREF ( K ) -XS ( I T ) ) *SIN(SAI ) 

X (K ) =XS (IT ) ♦ (XREF (< ) -XS ( IT ) ) *COS ( S A I ) ♦ ( YR£F ( K ) -YS ( I T ) ) »S IN ( SA I ) 
GO TO 21 

determine The trailing edge intersection 

between fixed and variable s»/eep whg sections 
IF ( ABS (KSAR (K ) -SAR (K-1 )) ,LT, (. 1/RAD) ) GO TO 32 
Y(K)=XREF(K*1)-X(K-1 ) -a (K) *YREF (<*1 )»A(K-1)«Y(K-1) 
Y(K)*Y(K)/(A(K-1)-A(K)) 

X(K)=A(K)*X(k- 1)-A(K-1) *XREF (K* 1 )*A(K-1)*A(K)* ( YREE (K» 1 ) -Y (K-l ) ) 
X(K)=X(K)/(A(K)-A(<-1) ) 

GO TO 33 
X(K)-XRET (K*i) 

Y (K ) =YREF (K* 1 ) 

K = K*1 

STORE reference PLANFORM COORDINATES ON INBOARD FIXED TRAll ING 
EDGE 

DO 39 I=K,N1 
X(I)=XREF(I) 

T(I)=YREF(I) 

SAR ( l-l ) =RSAR( l-l ) 

00 36 1 = 1, N 
XX ( I , IT) =X ( I ) 

YY(I,IT)=Y(I) 

MMCD(I*IT)=MC0(I,IT) 

TTWD( I , I T) =DIH( I , IT ) 

AS(I,IT)=A(I) 

XX(N1,IT)=X(n1) 

yy(ni,it)=y(n1) 

AN ( IT) =AAN (IT) 

CONTINUE 

LINE uR BREAKPOINTS AMONG PLANFORMS 

BOTSV(1)=80TSV(2)=0. 

WRITE (6*108) 

00 99 1T=1 * IPLAN 

NIT = AN ( I I) ♦ 1 

DO 93 ITT=1, IPLAN 

IF ( ITT.EO. IT) GO TO 93 

NITT = AN( ITT) *1 

DO 92 I=1*NITT 

JPSV=0 

DO 38 JP=1*NIT 

IF ( YY ( JP, IT) .EO.YY ( I , ITT) ) GO TO 92 

CONTINUE 

00 39 JP=1*NIT 

IF ( YY ( jP* IT) .LT.YY ( I, ITT) ) GO TO 90 

CONTINUE 

GO TO 92 

IF(JP.EO.I) GO TO 92 
JPSV=JP 

IND=NIT-(JPSV-1) 

DO 91 JP=1 * I NO 


GEO 181 
GEO 182 
geo 183 
geo 189 
GEO 185 

geo 186 

GEO 187 

geo 188 

GEO 189 
GEO 190 
GEO 191 
geo 192 
GEO 193 
GEO 199 
GEO 195 
GEO 196 
GEO 197 
GEO 198 
geo 199 

geo 200 

GEO 201 
GEO 202 
GEO 203 
GEO 209 
GEO 205 
GEO 206 
GEO 207 
GEO 208 
GEO 209 
GEO 210 
GEO 211 
GEO 212 
GEO 213 
GEO 219 
GEO 215 
GEO 216 
GEO 217 
GEO 218 
GEO 219 
GEO 220 
GEO 221 
GEO 222 
GEO 223 
GEO 229 
GEO 225 
GEO 226 
GEO 227 
GEO 228 
GEO 229 
GEO 230 
GEO 231 
GEO 232 
GEO 233 
GEO 239 
GEO 235 
GEO 236 
GEO 237 
GEO 238 
geo 238A 
GEO 239 
GEO 290 
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K2=NIT-JP*<? 



geo 

241 


K1=NIT-JP*1 



geo 

242 


XX(K2»IT)=XX(k1»IT) 



GEO 

243 


YY(K2.IT)=YY(K1«IT) 



geo 

244 


MMCD(K2» IT) =MMC0(K1* IT) 



geo 

245 


AS(K2*1T)=AS(KI.IT) 



GEO 

246 

41 

TTrtO<K2.IT)=TTWO(Kl.IT) 



GEO 

247 


YY(JPSV*IT)=YY(I»ITT) 



GEO 

248 


AS( JPSV»IT)=AS(JPSV-1»IT) 



GEO 

249 


TTWD( JPSV,IT)=TTWO<JPSV-l. 

IT) 


GEO 

250 


XX< JPSV» 1T)=(YY( JPSV»IT)-YY(JPSV-1*IT) )*AS( JPSV-l» 

IT)»XX(JPSV-1»ITGE0 

251 


1) 



GEO 

252 


MMCO( JPSV* IT) =MMCO( JPSV-1 * 

IT) 


GEO 

253 


AN(IT)=AN(IT)*1, 



geo 

254 


NIT=NIT*l 



GEO 

255 

42 

CONTINUE 



GEO 

256 

43 

CONTINUE 



geo 

257 

C 




GEO 

258 

c 

SEQUENCE WING COORDINATES 

FROM TIP TO ROOT 


GEO 

259 

c 




GEO 

260 


N1=AN(IT)*1. 



GEO 

261 


00 44 1=1. N1 



GEO 

262 

44 

Q(I)=YY(I.IT) 



GEO 

263 


DO 48 Jsl.Nl 



GEO 

264 


HIGH=1 . 



GEO 

265 


DO 45 I=1.N1 



geo 

?66 


IF ( (Q(D-HIGH) .GE.O.) GO 

TO 45 


geo 

267 


HIGH=Q ( I ) 



GEO 

268 


IH=I 



GEO 

269 

45 

CONTINUE 



GEO 

270 


IF (J.NE.l) GO TO 46 



GEO 

271 


80TSV ( IT) =HIGH 



GEO 

272 


KFX ( IT) =IH 



geo 

273 

46 

Q(IH)=1. 



GEO 

274 


SPY(J.IT)=HIGH 



geo 

275 


IF (IH.GT.KFX(IT) ) GO TO 47 


geo 

276 


IYL(J.IT)=1 



geo 

277 


IYT(J.IT)=0, 



GEO 

278 


60 TO 48 



GEO 

279 

47 

lYL ( J. IT) =0 



geo 

280 


IYT(J»IT)=1 



GEO 

281 

48 

CONTINUE 



geo 

282 

49 

CONTINUE 



geo 

2H3 

C 




geo 

284 

C 

SELECT MAXIMUM B/2 AS THE 

WING semispan. if 

both first And geo 

285 

c 

SECOND PLANFORMS HAVE SAME SEMlS=AN THEN THE 

SECOND PLANFORM IS GEO 

286 

c 

taken to be THE WING, 



geo 

287 

c 




GEO 

288 


KB0T=1 



GEO 

289 


IF (BOTSV(l) ,GE.B0TSV(2) ) 

KB0T=2 


GEO 

?90 


B0T=8OTSV (KBOT) 



geo 

291 

c 




GEO 

292 

c 

COMPUTE nominal HORSESHOE 

VORTEX width along 

MING 

SURFACE GEO 

293 

c 




GEO 

294 


TSPAN=0 



GEO 

295 


ISAVE=KFX (K80T)-1 



geo 

296 


I=KFX(KB0T)-2 



GEO 

297 

50 

IF (I.EO.O) GO TO 51 



GEO 

298 


IF (TTWO( I.KBOT) .EQ.TTWDdSAVE.KSOT) ) GO TO 52 

GEO 

299 

51 

CTWD=C0S( ATAN(TTW0( ISAVE. 

KBOT) ) ) 


GEO 

300 
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52 

53 
C 


54 

55 


56 

57 
C 
C 

c 


5B 

59 

C 

c 

c 


60 

51 


52 


53 

C 

C 


TL6TH= ( YY ( ISAVE*1«<B0T)-YY( 1*1 ,K80T) ) /CTWD 
TSPAN=TSPAN*TLGTh 
IF (I.EQ.O) GO TO 53 
ISAVE=I 
1 = 1-1 
GO TO 50 
VI=TSPAN/VIC 
VSTOL=VI/2 

TSPANA=0. 

KBIT=2 

IF (IPLAiN.EQ.l) GO TO 57 

IF (K80T.E0.2) KBIT=1 

1SAVEA=KFX(KBIT)-1 

IA=KFX(KH1T)-2 

IF ( lA.EO.O) GO TO 55 

IF (TTlYDI IA«KHIT) .EO.TTWOI ISAVFA.kBIT) ) GO TO 56 
CTrtOA=COS(ATAN<TTWO( ISAVEA*KBIT))) 

TLGTHA=(YY( ISAVEA*1»KBIT)-YY( IA*1.KBIT) )/CTwOA 

T3PANA = 7SPAf'jA^TLGT HA 

IF ( lA.EO.O) GO TO 57 

ISAVEA=IA 

IA=IA-1 

GO TO 54 

CONTINUE 

Eliminate planform breakpoints which are within (b/2)/20oo units 
laterallt 

DO 59 IT=1*IPLAN 
N=AN ( IT ) 

N1=N*1 
DO 59 J=1,N 

AA = ABS(SPY ( J, IT) -SOY ( J*1 ♦ IT) ) 

IF ( AA.EO.O . .09. AA.GT.ABS (TSPAn/2000. ) ) GO TO 59 
IF (AA.gT.YTOL) write (6«111) spy { J» 1 . it ) .spy ( J, IT) 

DO 58 I=1.N1 

IF (YY(I.IT).NE.SPY(J»1,IT)) go to 58 
YY(I.IT)=SPY(J,IT) 

CONTINUE 

SPY(J*1.IT)=SPY(J»IT) 

CONTINUE 

COMPUTE L COORDINATES 

DO 63 IT=1.IPLAN 
UM=N1=AN( IT) *1. 

DO 60 JZ=1 *N1 
ZZ(jZ.IT)=RTCDHT(IT) 

JZ=1 

UZ=JZ*1 

IF ( JZ.GT.KFXI IT) ) GO TQ 62 

ZZ{ JZ. IT ) =ZZ ( jZ-lf IT) ♦ (YY ( JZ. IT)-YY ( JZ-1 . IT) ) " TT WO ( JZ- 1 ♦ I T ) 

GO TO 61 
JM=JM-1 

IF ( JM.EO.KFX(IT) ) GO TO 63 

ZZ(JM.IT)=ZZ(JM*1*IT)*(YY(JM.IT)-YY(JM*1.IT))«TTW0(JM,IT) 

GO TO 52 
CONTINUE 


GEO 301 
GEO 302 
GEO 303 
GEO 304 
GEO 305 
GEO 306 
GEO 307 
GEO 308 
GEO 309 
GEO 310 
GEO 311 
GEO 312 
GEO 313 
GEO 314 
GEO 315 
GEO 316 
GEO 317 
GEO 318 
GEO 319 
GEO 320 
GEO 3Pl 
GEO 322 
geo 323 
GEO 324 
GEO 325 
GEO 326 
geo 327 

geo 328 

GEO 329 
GEO 330 
GEO 331 
GEO 332 
geo 333 
GEO 334 
GEO 335 
GEO 336 
GEO 337 
GEO 338 
GEO 339 
geo 340 
GEO 341 
GEO 342 
GEO 343 
GEO 344 
GEO 345 
GEO 346 
GEO 347 
GEO 348 
geo 349 
GEO 350 
GEO 351 
GEO 352 
GEO 35-3 
GEO 354 
GEO 355 
GEO 356 
GEO 357 
GEO 3S8 


WRITE planform PERIMETER POINTS 


GEO 359 

actually used in THE COMPUTATIONS GEO 360 
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c 


65 

C 

c 

c 


c 

66 


C 

c 


c 

67 
C 

c 

c 

68 


69 


70 


71 


C 

72 


GEO 

WRITE (6*100) 

DO 65 IT = 1*IPL^N 

n=an(IT) geo 

N1=N*1 geo 

IF (IT,e0.2) WRITE (6*110) GEO 

DO 64 KK=1*N geo 

TOUT=ATAN(TTWO(KK* IT) )«RA0 GEO 

A0UT=ATAM(AS(KK*IT) )*RAD GEO 

IF (AS(KK»IT) .EQ.AZY) AOUT=90. GEO 

WRITE (6*101) KKfXX(KK*lT)*YY(KK*lT)*ZZ(KK*lT) * AOUT *TOUT *MMC0 ( KK» I GEO 
IT) GEO 

CONTINUE geo 

WRITE (6*101) N1*XX(N1*IT)*YY(n1*IT)*ZZ(N1*IT) GEO 

CONTINUE GEO 

GEO 

PART ONE - SECTION THREE - LAY OUT YAWED HORSESHOE VORTICES GEO 

GEO 
GEO 

geo 


STRUE=0. 

NSSWSV ( 1 ) =NSSWSV(2) =MSV ( I ) =MSV (2) =0 


DO 74 IT=1*IPLAN 
N1=AN( IT) ♦! . 

1 = 0 
J=1 

YIN=B0TSV (IT) 

ILE=ITE=KFX( IT) 

determine SPANWISE borders of horseshoe VORTICES 

IXL=IXT=0 

I=I^1 

IF (YIN.GE. (SPY(J* IT) ♦VSTOL) ) GO TO 67 

BORDER IS WITHIN VORTEX SPACING TOLERANCE (VSTOL) 

therefore USE THE NEXT BREAKPOINT INBOARD FOR THE 

VB0R0(I)=YIN 

GO TO 70 

USE NOMINAL VORTEX SPACING TO DETERMINE THE BORDER 
VBORO( I ) =SPY ( J* IT) 

COMPUTE SUBSCRIPTS ILE AND ITE TO INDICATE WHICH 
BREAKPOINTS ARE ADJACENT AND WHETHER THEY ARE ON 
EDGE OR THF trailing EDGE 
IF (J.GE.Nl) GO TO 69 

IF (SPY(J*IT) .NE.S=>Y(J^1*IT) ) 30 TO 69 
IXL=IXL^1YL( J*IT) 

IXT=IXT^IYT ( J*IT) 

J=J^1 
60 TO 60 
YIN=SPY( J*IT) 

IXL=IXL*IYL(J*IT) 

IXT=IXT^IYT(J*IT) 

J=J^1 

CPHI=C0S(ATAN(TTWD( IlE* IT) ) ) 

IPHI=ILE-IXL 
IF (J.GE.Nl) IPHI=1 

YIN=YIN-VI»C0S (ATAN (TTWD( IPHI* IT) ) ) 

ir (I.NE.l) GO TO 72 

ILE=ILE-IXL 

ITE=ITE^IXT 

60 TO 66 

COMPUTE COORDINATES FOR CHOROWISE ROW OF HORSESHOE 
YO=(VBORO( I-l ) ♦VBORO( I ) )/2. 

HW= ( VBORD ( I ) -VBORD ( I - 1 ) ) /2. 


OF BREAKPOINT 

border 


GEO 
GEO 

geo 
geo 
geo 
geo 

GEO 

geo 
geo 

GEO 
GEO 
GEO 
GEO 

geo 

GEO 
GEO 
GEO 

THE WING leading GEO 

geo 

GEO 

geo 

geo 

geo 

geo 

geo 

geo 

geo 

geo 

GEO 

geo 

geo 

geo 

geo 

GEO 

geo 

geo 

GEO 

VORTICES GEO 

GEO 
GEO 


361 

362 

363 

364 

365 

366 
367 
360 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 

380 
381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 
3^4 

395 

396 

397 

398 

399 

400 

401 
A02 

403 

404 

405 

406 

407 

408 

409 

410 

411 

412 

413 

414 

415 

416 

417 

418 

419 
4?0 
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73 

C 

C 

C 

It* 

C 

c 

c 


7S 


1M1 = I-1 ♦NSSWSVn ) 

i!H(IMl)=Z2(lLE.IT)*(YQ-yY(ILE,lT) ) *TTWD ( I LE « I T ) 

PHI ( IMl ) =TTWO( ILE«IT) 
sswwA(ih1)=as(ile«iT) 

XlE=XX(ILE.IT)*AS(ILE.IT)*(YO-yY(ILE*IT) ) 

XET=XX( ITE» IT) ♦AS(ITE*IT)*(YO-yY( ITE* IT) ) 
XLOCAL=(XLE-XET)/TBLSCtt( IMl) 

COMPUTE WING AREA PROJECTED TO The X - Y PLANE 

STRUE=STRUE*XL0CAL*T8LSCW‘( IMl ) *(HW*2. ) »2. 

NSCW=T8LSCW ( IMl ) 

00 73 JCW=1,NSCW 
AJCW=JCW-1 

XLEL=XLE-AJCW*XLOCAL 
NTS=JCW*MSV ( 1 ) «MSV (2) 

PN(NTS)=XLEL-.25«XlOCAL 

PV(NTS)=XLEL-.75*XL0CAL 

PS I (NTS) = ( ( XLE-PN (NTS) ) * AS ( I TE ♦ I T ) ♦ ( PN ( NTS ) -XET ) * AS ( ILE ♦ 
IXET) *CPHl 
S(NTS) =HW/CPHI 
0(NTS)=YQ 
CONTINUE 

MSV (IT) =MSV ( IT) *NSCW 


TEST TO DETERMINE WHEN WING ROOT 
IF ( V80R0( I ) ,LT.YREG( 1 * IT) ) GO TO 


a-1 


IS reached 

71 


7b 


NSSWSV( IT); 

CONTINUE 
M=MSV ( 1 ) *MSV (2) 

COMPUTE aspect ratio AND AVERAGE CHORD 
bOT=-0OT 

AR=A.*B0T»B0T/SREF 
ARTRUE=4.*BOT*flOT/STRUE 
CAVE=STRUE/ (2.*B0T) 

BET A= ( 1 ,-MACH^MACH) **.5 
WRITE (6*11A) M 

WRITE (6*11S) ( IT*MSV ( IT) »NSSWSV ( IT ) ♦ IT=1 . IPLAN) 
IF (SCW.NE.O.) WRITE (6fll2) SCW 

IF (SCW.EQ.O.) WRITE (6.113) ( TBLSCW ( I ) . 1 = 1 , NSTA ) 

apply PRANOTl-GLAUERT CORRECTION 
DO 75 NV=1.M 

PSI (NV)=ATAN(PSI (NV)/8ETA) 

PN(NV)=PN(NV)/BETA 
PV(NV)=PV(NV)/HETA 
NSSW=NSSWSV ( 1 ) .NSSWSV (2) 

JN = 0 

00 77 JSSw=l,NSSW 
CH0RD(JSSW)=0. 

NSCW=TBlSCw( JSSW) 

DO 76 JSCW=1,NSCW 
JN=JN*1 

ChORD(JSSW) = CHORD ( JSSW )-2.*(?V(JN)-PN(JN) )*HETA 
CONTINUE 


GEO 421 
GEO 422 
GEO 423 
GEO 424 
GEO 425 
GEO 426 
GEO 427 
GEO 42B 
GEO 429 
GEO 430 
GEO 431 
GEO 432 
GEO 433 
GEO 434 
GEO 435 
GEO 436 
GEO 437 
GEO 438 
GEO 439 
IT) )/(XLE-GE0 440 
GEO 441 
GEO 442 
GEO 443 
GEO 444 
geo 445 
GEO 446 
GEO 447 
GEO 448 
GEO 449 
geo 450 
GEO 451 
GEO 452 
GEO 453 
GEO 454 
GEO 455 
GEO 456 
GEO 457 
GEO 458 
GEO 459 
GEO 460 
GEO 461 
GEO 462 
GEO 463 
GEO 464 
geo 465 
geo 466 
GEO 467 
GEO 468 
GEO 469 
GEO 470 
GEO 471 
GEO 472 
GEO 473 
GEO 474 
GEO 475 
GEO 476 
GEO 477 
geo 478 
GEO 479 
geo 480 
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77 

XTE ( JSSW) = (PV ( JN) ♦ (PV < JN)-PN( JN) ) /2, ) *BETA 

geo 981 


PHISUM=0. 

geo 982 


00 78 IKY=1,NSSW 

geo 983 


PHISUM=PHISUM+PHI (IKY) 

geo 989 

78 

CONTINUE 

geo 985 


IFLA6=I 

geo 986 


IF (IPLAN.EO.I.ANO.PHISUM.NE.O.) IFLA6=2 

geo 987 


IF ( IPLAN.E(3.2. ANO.PhISUM.NE.O. ) GO TO 79 

geo 988 


GO TO 83 

geo 989 

79 

00 81 IP=1.IPLAN 

geo 990 


IA=1»(IP-1)*NSSWSV(1) 

geo 991 


IB=NSSWSV ( 1 ) ♦ ( lP-1 ) “NSSWSV (2) 

geo 992 


IC=1-(IP-2)*NSSWSV(1) 

geo 993 


IO=NSSWSV ( 1 ) - ( IP-2) *NSSWSV (2) 

geo 499 


00 80 IU=IA»18 

GEO 995 


00 80 IZ=IC*IO 

geo 996 


IF (ZH(IU) .EQ.ZH(IZ)) go to 82 

geo 997 

80 

CONTINUE 

geo 998 

81 

CONTINUE 

geo 999 


IFLAG=2 

geo 500 


GO TO 83 

geo 501 

82 

IFLAG=3 

geo 502 

83 

CONTINUE 

geo 503 


REAO (5»122) XCFWtXCFT 

geo 509 


IF (N.GT.AOO) GO TO 86 

geo 505 


NSW=NSSWSV(1) ♦NSS»(SV(2) 

geo 506 


IF (NSW.GT.50) GO TO 85 

GEO 507 


ITSV=0 

GEO 508 


00 89 IT=1.IPLAN 

geo 509 


IF (AN(IT) .LE.25.) GO TO 89 

GEO 510 


WRITE (6*118) IT,AN(1T) 

GEO 511 


ITSV=1 

GEO 512 

84 

CONTINUE 

geo 513 


IF (ITSV.GT.O) go to 91 

GEO 519 


GO TO 87 

geo 515 

3S 

WRITE (6*117) NSW 

GEO 516 


GO TO 91 

geo 517 

86 

WRITE (6*116) M 

GEO 518 


GO TO 91 

GEO 519 

87 

REWIND 25 

GEO 520 


WRITE (25) e0T*M*8ETA*PTEST*0TEST*TBLSCW*Q*PN.PV*S*PSI*PHl 

,ZH*NSSWGE0 521 


1 * TWIST *CREF.SRFF* CAVE* CLOES*STRUF*AR *ARTRuE*RTCDHT* CONFIG* 

NSSwSVtMGEO 522 


2SV*KB0T*PLAN* IPUAN*MACH*SSWWA*CH0R0*XTE«KRIT*TSPAN*TSPANA* 

XCFW*XCFGE0 523 


3T , IFL AG* YREG (1*1)* YREG (1*2) 

geo 529 


END file 25 

geo 525 


GO TO (88*89*90) * IFLAG 

geo 526 

88 

WRITE (6*119) 

geo 527 


WRITE (50*123) 

GEO 528 


GO TO 92 

GEO 529 

89 

WRITE (6*120) 

geo 530 


WRITE (50*129) 

geo 531 


GO TO 92 

geo 532 

90 

WRITE (6*121) 

geo 533 


WRITE (50*125) 

geo 539 


GO TO 92 

geo 535 

91 

T0TAL=T0TAL-1 . 

GEO 536 


WRITE (50*126) 

geo 537 

92 

CONTINUE 

geo 538 


END file 50 

geo 539 


STOP 

GEO 590 
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93 

ICODEOF=l 

GEO 

541 


WRITE (6.103) CONFIG 

GEO 

542 


GO TO 91 

GEO 

543 

94 

ICODEOF=2 

GEO 

544 


WRITE (6*104) k.IT 

GEO 

545 


GO TO 91 

GEO 

546 

9b 

ICODEOF=3 

GEO 

547 


WRITE (6*107) PTEST, OTEST 

GEO 

540 


GO TO 91 

GEO 

549 

C 


GEO 

550 

C 


GEO 

551 

96 

FORMAT ( 1H1//63X, 13HGE0METRY DATA) 

GEO 

552 

97 

FORMAT (///45X.A10.22HREFERENCE pl ANFORM HAS.I3.7H CUR VES// 1 2X , 1 9HGE0 

553 


IROOT CHORD HEIGHT = ,F 1 2 . 5.4X ,?9HV AR I ABLE SWEEP PIVOT 

position. 4X.6GE0 

554 


2HX(S) =.F12.5.5X.6HY (S) = ,F 1 2 . 5//46X * 40HBRE AK POINTS 

FOR The refergeo 

555 


3ENCE PLANFORM /) 

geo 

556 

90 

FORMAT (8F10.4) 

geo 

557 

99 

FORMAT (1H1//47X.17HCONFIGURATION NO..F8.0/) 

GEO 

558 


100 

101 

10 ^ 

103 

104 

105 

106 

107 

lOB 

109 

110 
111 

112 

113 

114 

115 

116 

117 

IIB 

119 

120 
121 
122 

123 

124 

125 

126 


FORMAT (22X .5HP0INT«6X« IHX, 1 lx. IHY. 1 lx. 1H2. 1 OX . 5MSWEEP* 7X . 8H0I HtOHGEO 559 
iAL.4X.4HM0vt/6BX.5HAMGLF.flX,5HANGLE.6X.4HC0DE/) GEO 560 

FORMAT (20X.I5.3F12.5.2F14.5.I6) GEO 56l 

FORMAT (/40X.5HCUW0E. I3.9H IS SWEPT . F 1 2 . 5 . 20H DEGREES ON PLANFORM.GEO 562 
1^3) f-jEO 563 

FORMAT ( 1h1///41X.43hEN0 OF FILE ENCOUNTERED AFTER CONF I GURAT 1 ON.FGEO 564 

GEO 565 

FORMAT ( 1H1///18X.45HTHE FIRST VARIABLE SWEEP CURVE SPECIFIED (K =gEO 566 
1.I3.44H ) DOES NOT HAVE AN M CODE OF 2 FOR PLANFORM.I4) GEO 567 

FORMAT(5F5.1.2F10.4) qEq 568 

FORMAT (26X.I5.2F12,5.2F16,5.4X.T4) GEO 569 

FORMAT ( 1H1///1X.38HERROR - PROGRAM CANNOT PROCESS PTEST =.F5,1.12GE0 570 
IH AND OTEST =.F5.1) GEO 571 

FORMAT (//48X.35H8REAK POINTS FOR THIS CONFIGURATION//) GEO 572 

FORMAT ( 28X.5HPOINT . 6X . 1 hX . 1 1 x . 1 HY . 1 1 X . bHSwEEP . 1 OX . BHDIHEOR AL . 7X ,4 gEO 5 73 
1HMOVE/3BX. 3HREF »9X.3hREF . 10X.5H ANGLE. 1 1 X ♦ 5H ANGLE » 9X . 4HC00E/ ) GEO 574 

FORMAT (/52X,28HSEC0ND PLANFORM BREAK POINTS/) GEO 575 

FORMAT (////25X.34HTHE BREAKPOINT LOCATED SPANWISE AT . F 1 1 .5 . 3X . 20HGEO 576 
IHAS BEEN ADJUSTED T0.F9.5////) GEO 577 

FORMAT (/43X.F5,0.41h HORSESHOE VORTICES IN EACH CHORDWISE ROW) GEO 578 
FORMAT (/23X.9BHTA3LE OF HORSESHOE VORTICES IN EACH CHORDWISE ROW GEO 579 
KFROM TIP TO ROOT BEGINNING WITH FIRST PL ANFORM ) //25F5 . 0/25F5 . O ) GEO 5B0 
FORMAT (///33XI5.62H HORSESHOE VORTICES USFD ON THE LEFT HALF OF TGEO 581 
IHE C0NFIGURATION//50X.36HPLANFORM TOTAL SPANWISE/) GEO 582 

FORMAT (p 2X. 14, lOX.n.llX. 14) GEO 5B3 

FORMAT ( iHl/ZlOX. I6.93HH0RSESH0E VORTICES LAIDOUT. THIS IS MORE THGEO 584 
IAN The 400 MAXIMUM. THIS CONFIGURATION IS ABORTED.) GEO 585 

FORMAT (1H1//10X.I6.101H ROWS OF HORSESHOE VORTICES LAIDOUT. THIS GEO 586 
IIS MORE Than the 50 maximum, this CONFIGURATION IS ABORTED.) GEO 587 

FORMAT ( 1h1//10X»8hPlANF0RM, I6.4H HAS.I6.74H BREAKPOINTS. THE MAxIGEO 588 


IMUM DIMENSIONED IS 25. THE CONFIGURATION IS ABORTED.) 
format (///20X,28H'4InIMUM FIELD LENGTH = 51000) 
format (///20X.28HMINIMUM FIELD LENGTH = 63000) 
format (///20X.29HMIN1MUM FIELD LENGTH = 112000) 
format (6F10.4) 

FORMAT (J5h*C0MP1LE VLMCDR AGS. VL MCC I R 1 , VLMC20C ) 

FORMAT (35h*C0MPILE VLMCDRAGS , VL MCC 1 R2 . VLMCZOC ) 
format (35H*C0MPILE VLMCDR AGS . VL MCC I R3 . VLMCZOC ) 

FORMAT (18H«C0MPILE VLMCDUMMY) 

END 


geo 589 


GEO 

GEO 

GEO 

GEO 

GEO 

GEO 

GEO 

GEO 

GEO 


590 

591 

592 

593 

594 

595 

596 

597 
598- 
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»OECK VLMCD«AGS 

OVERLAY twlN6TL*0*0) 060 

PROGRAM XIN6AL (OUTPUT. TAPE6=OlJTPOT»TAPE10»TAPE20«TAPE2!>) 060 

COMMON /ALL/ BOT *M* BET A , PTEST ♦ QTEST ♦ TBlSCw ( 50 ) » 0 (400 ) » PN (400 ) .PV (4060 
100) *S(400) .PS I (400) .PHI (50) .ZH(50) .NSSi( 060 

COMMON /ONEThRE/ Ti^IST ( 2) .CREF.SREE.CAVE.CLDES.STRUE. AR.AHTRUE.RTC060 
lOHT (2) .C0 NFIG.NSSi«(SV (2) .MSV (2) .KBOT .PLAN. I PL AN. MACH. SS<(W A (50) .XCF«OGO 
2.XCFT.YREG( 1 .2) 060 

COMMON /TOTHRE/ CIR(400) 060 

COMMON /CCRRDO/ CHORD (50 ) .XTE ( 50 ) .KBIT. TSPAN . TSPANA 060 


real MACH 


VORTEX LATTICE AERODYNAMIC COMPUTATION 
REWIND 2i> 

READ (25) BOT.m.BETA.PTEST.QTEST.TBLSCw.Oi 


3 .IFLAG.YREG(l.l) .YREG(1.2) 


WIN6TL=6LWINGTL 

RECALL=6HRECALL 


CALL overlay (WINGTL. 1.0. recall) 


060 


CALL 

STOP 

END 


overlay (WINGTL. 2.0. RECALL) 


1 

2 

3- 

4 

5 

6 

7 

8 
9 


DGO 

10 

OGO 

11 

OGO 

12 

DGO 

13 

DGO 

14 

DGO 

15 

OGO 

16 

DGO 

17 

DGO 

18 

OGO 

19 

060 

20 

OGO 

21 

DGO 

22 

060 

23 

OGO 

24 

060 

25 

DGO 

26 

DGO 

27 

DGO 

28 

DGO 

29 

DGO 

30 

DGO 

31- 
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SUBROUTINE FTLUP ( X ♦ Y * ♦ N« VAR I , V ARO ) 


TLU 

1 

***D0CUMENT date 09-12-69 SUBROUTINE REVISED 

07-07-69 

2 

MODIFICATION OF LIBRARY INTERPOLATION 

SUBROUTINE FTLUP TLU 

3 

DIMENSION VARKD* VARD(l)* V(3)» YY(2) 


TLU 

4 

DIMENSION 11(43) 


TLU 

5 



TLU 

6 

initialize all interval pointers to -1.0 

FOR 

monotonicity checktlu 

7 

DATA (II ( J) .J=I.43)/43*-l/ 


TLU 

8 

MA=IABS(M) 


TLU 

9 



TLU 

10 

ASSIGN INTERVAL POINTER FOR GIVEN 

VARI TABLE TLU 

11 

THE SAME POINTER RILL BE USED ON A GIVEN 

VARI 

TABLE EVERY TIME TLU 

12 

LI=MOD(LOCF (VARI (1 ) ) ,43) ♦! 


TLU 

13 

I = II (LI ) 


TLU 

14 

IF (I.GE.O) GO TO 6 


TLU 

15 

IF (N.LT.2) GO TO 6 


TLU 

16 



TLU 

17 

MONOTONICITY CHECK 


TLU 

18 

IF < VART ( 2) -VAQI ( 1 ) ) ?.2.4 


TLU 

19 

ERROR IN monotonicity 


TLU 

20 

K=L0CF (VARI ( 1 ) ) 


TLU 

21 

PRINT 17» J,K. (VARI ( J) t J=1.N) , (VAR0( J) * J= 

1 ,N) 

TLU 

22 

STOP 


TLU 

23 

MONOTONiC DECREASING 


TLU 

24 

DO 3 J=2»N 


TLU 

25 

IF (VARI ( J)-VARI(J-l) ) 3«1«1 


TLU 

26 

CONTINUE 


TLU 

27 

GO TO 6 


TLU 

28 

monotonic increasing 


TLU 

29 

DO 5 J=2»N 


TLU 

30 

IF (VARI ( J)-VARI (J-1) ) 1,1.5 


TLU 

31 

CONTINUE 


TLU 

32 



TLU 

33 

interpolation 


TLU 

34 

IF (I.LE.O) 1=1 


TLU 

35 

IF (I.GE.N) I=N-1 


TLU 

36 

IF (N.LE.1) go to 7 


TLU 

37 

IF (MA.nE.O) go to 8 


TLU 

38 

ZERO ORDER 


TLU 

39 

Y=VARD( 1 ) 


TLU 

40 

GO TO 16 


TLU 

41 

LOCATE I INTERVAL < X ( I ) . LE ,X .L T . X ( I ♦ 1 ) ) 


TLU 

42 

IF ( ( VARI ( I) -X)*(VARI ( 1,1 )-X) ) 11,11,9 


TLU 

43 

IN GIVES DIRECTION FOR SEARCH OF INTERVALS 

TLU 

44 

IN = SIGN(1.0, (VARI (1*1) -VARI (1) )*(X-VARI (I) ) ) 

TLU 

45 

IF X OUTSIDE ENDPOINTS, EXTRAPOLATE FROM 

END INTERVAL TLU 

46 

IF ( (I*IN) ,LE.O) GO TO 11 


TLU 

47 

IF ( ( I»IN) .GE.N) GO TO 11 


TLU 

48 

I=1*IN 


TLU 

49 

IF ( ( VARI ( I )-X)*(VARl ( 1*1 )-X) ) 11,11,10 


TLU 

50 

IF (MA.EO.2) go to 12 


TLU 

51 



TLU 

52 

FIRST ORDER 


TLU 

53 

Y=(VARD( I)*(VARI (I,1)-X)-VARD(I»1 )*(VARI ( I ) -X ) ) / ( V AR I < I ♦ 1 ) - V AR K I ) TLU 

54 

1 ) 


TLU 

55 

GO TO 16 


TLU 

56 



TLU 

57 

SECOND O^OEW 


TLU 

58 

IF (N.E0.2) GO TO 1 


TLU 

59 

IF (I,EO. (N-1) ) 60 TO U 


TLU 

60 
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ir (I.EO.l) 60 TO 13 

TUU 

61 

c 

PICK thi«o point 

TLU 

62 


SK=VARI ( !♦! ) -VARI ( I ) 

TLU 

63 


IF ( <SK* (X-VARI ( I-l ) ) ) .LT. (SK*(VARI (I*2)-X) ) ) 60 TO 14 

TLU 

64 

13 

L=I 

TLU 

65 


60 TO 15 

TLU 

66 

14 

L=I-1 

TLU 

67 

lb 

V(1)=VARI (U-X 

TLU 

68 


V(2)=VARI (L*l )-X 

TLU 

69 


V(3)=VAR1 (L*2)-X 

TLU 

70 


rv(i ) = <varo(l)*v(2)-varo(L*i )*v(i ) )/(vari (l*i )"Vari (d > 

TLU 

71 


TY<2) = ( VAR0(L* 1 ) *v (3>-VAR0(U*2)*V(2) )/(VARi (L*2)-VARKL*1)) 

TLU 

72 


Y=(YV(l )*V(3)-YY(2)*V< 1) )/(VARI (L *2) -VARI (L) ) 

TLU 

73 

16 

II(LI)=I 

TLU 

74 


RETURN 

TLU 

75 

C 


TLU 

76 

C 


TLU 

77 

17 

FORMAT (1H1.50H TABLE BELOW OUT OF ORDER FOR FTLUP AT POSITION 

tTLU 

78 


1I5*/31H X table is STORED IN LOCATION ♦06*// (8615,8) ) 

TLU 

79 


END 

TLU 

80- 


121 



APPENDIX E 


c 

c 

c 


c 

c 

c 

1 


2 

c 

c 


3 

4 
3 

6 

7 

B 


9 
C 
C 
C 

10 


11 

12 


13 

14 

C 

C 

c 

15 

16 


17 


SUBROUTINE SIMF.O ( A, N»B » H t OETEBM . IP I VOT * NM AX . I SCALE ) 

SOLUTION OF simultaneous LINEAR EQUATIONS 

*** document date 08-01-68 SUBROUTINE REVISED 08-01-68 

DIMENSION IPlVOT(N)t A(NMAX*N), B(NMAX,M) 

equivalence (IR0W»JR0W)» ( ICDLJM» JCOLUM) » ( AMAX»T«SWAP> 

initialisation 

iscale=o 

Rl=10.0*»100 
R2=1.0/Rl 
DETERM=1.0 
DO 2 J=1«N 
IPIVOT ( J)=0 
DO 38 I=1.N 

SEARCH FOR PIVOT ELEMENT 

AMAX=0.0 
DO 7 J=liN 

IF ( IPIVOT U)-l ) 3.7»3 
DO 6 K=1 *N 

IF ( IPIVOT (K)-l ) 4.6.39 

IF (ABS(AMAX)-ABS(A( J.K) ) ) 5.6.6 

IROW=J 

ICOLUM=K 

AMAX=A(J.K) 

CONTINUE 

CONTINUE 

IF (AMAX) 9.8.9 

0ETERM=0.0 

ISCALE=0 

GO TO 39 

IPIVOT ( I COLUM)= IPIVOT ( ICOLUM) *1 

interchange ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

IF ( IROW-ICOLUM) 10.14.10 
OETERMs-DETERM 
DO 11 L=1«N 
SWAP=A ( IROW.L) 

A(IROW.L)=A(ICOLUM.L) 

A ( ICOLUM. L ) =SWAP 
IF (M) 14.14.12 

DO 13 L=1.M 
SWAP = 8 ( HOW.L) 

B(IR0W.L)=B<1C0LUM.L) 

B( ICOLUM.L)=SWAP 
Pi V0T= A ( ICOLUM. ICOLUM) 

IF (PIVOT) 15.8.15 

scale the determinant 

PI V0TI=PIV0T 

IF (ABS(OETERM)-Rl) 18.16.16 

DETERM=D£TERM/R1 

ISCALE=ISCALE»1 

IF (ABS(OETERM)-Rl) 21.17.17 

DETERMsDETERM/RI 


SEO 

1 

SEO 

2 

* 5 £ Q 

3 

SEQ 

4 

SEQ 

5 

SEO 

6 

SEQ 

7 

SEQ 

8 

SEQ 

9 

SEQ 

10 

SEQ 

11 

SEQ 

12 

SEO 

13 

SEQ 

14 

SEQ 

15 

SEQ 

16 

SEQ 

17 

SEQ 

10 

SEQ 

1 o 

SEQ 

20 

SEQ 

21 

SEQ 

22 

SEQ 

23 

SEQ 

24 

SEQ 

25 

SEQ 

26 

SEQ 

27 

SEQ 

28 

SEO 

29 

SEQ 

30 

SEQ 

31 

SEQ 

32 

SEQ 

33 

SEQ 

34 

SEQ 

35 

SEQ 

36 

SEQ 

37 

SEQ 

38 

SEO 

39 

SEQ 

40 

SEQ 

<41 

SEO 

42 

SEQ 

43 

SEQ 

44 

SEQ 

45 

SEO 

46 

SEO 

47 

SEQ 

48 

SEQ 

49 

SEQ 

50 

SEO 

51 

SEQ 

S2 

SEQ 

S3 

SEQ 

54 

SEQ 

55 

SEQ 

56 

SEO 

57 

SEQ 

58 

SEQ 

59 

SEQ 

60 
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ISCALE=ISCALE*1 

SEO 

61 


GO TO 21 

SEO 

62 

18 

IF (ABS(0ETEPM)-R2) 19«19»21 

SEO 

63 

19 

DETERM=D£TERM*91 

SEO 

64 


ISCALE=ISCALE-1 

SEO 

65 


IF (ABS(0ETERM)-R2) 20*20*21 

SEO 

66 

20 

oeterm=deterh*ri 

SEO 

67 


ISCALE=ISCALE-1 

SEO 

68 

21 

IF (ABS(PIVOTI )-Rl ) 29*22*22 

SEO 

69 

22 

PIVOTI=PIVOTI/R1 

SEO 

70 


ISCALE=ISCAUE*1 

SEO 

71 


IF (ABS(PIVOTI)-Rl) 27*23*23 

SEO 

72 

23 

PIV0T1=PIV0TI/R1 

SEO 

73 


ISCALE=ISCALF*1 

SEO 

74 


GO TO 27 

SEO 

75 

24 

IF (ABS<PIV0TI)-R2) 25*25*27 

SEO 

76 

25 

PIVOTI=PIVOTI*R1 

SEO 

77 


ISCALE*ISCAUE-1 

SEO 

78 


IF (ABS(PIVOTI )-R2) 26*26*27 

SEO 

79 

26 

PIV0TI=PIV0TI*R1 

SEO 

80 


ISCALE=ISCALE-1 

SEO 

81 

27 

determ=oeterm*pivoti 

SEO 

82 

C 


SEO 

83 

c 

DIVIDE PIVOT ROW BY PIVOT ELEMENT 

SEO 

84 

c 


SEQ 

85 


DO 29 L=1*N 

SEO 

86 


IF (IPIVOT(L)-l) 28*29*39 

SEO 

87 

28 

A( IC0LUM*L) =A( ICOLJM*L)/PIVOT 

SEO 

88 

29 

CONTINUE 

SEQ 

89 


IF (M) 32*32*30 

SEQ 

90 

30 

DO 31 L*1*M 

SEQ 

91 

31 

B(ICOLUM*U=B( IC0LJM*L)/PIV0T 

SEO 

92 

C 


SEO 

93 

c 

REDUCE NON-PIVOT ROWS 

SEO 

94 

c 


SEQ 

95 

32 

DO 38 L1=1*N 

SEO 

96 


IF (Ll-ICOLUM) 33*38*33 

SEO 

97 

33 

T=A(L1*IC0LUM) 

SEQ 

98 


DO 35 L=1.N 

SEQ 

99 


IF (IPIVOT(L)-I) 39*35*39 

SEQ 

100 

34 

A(L1*U=A{L1*L)-A(IC0LUM*L)*T 

SEO 

101 

35 

CONTINUE 

SEQ 

102 


IF (M) 36*38*36 

SEO 

103 

36 

DO 37 L=1*M 

SEQ 

104 

37 

B(L1*L)=S(L1*L)-B( ICOLUM*L)*T 

SEO 

105 

38 

CONTINUE 

SEO 

106 

39 

RETURN 

SEQ 

107 


END 

SEQ 

108 
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SUBROUTINE DRaoSUB ( R . A» Y « Z ♦ S . I S« JS . WN< ) 

REAL IS.JS 

ZP=Z*S«SIN(A) 

YP=Y*S»C0S(A) 

ZM=Z-S»SIN(A) 

YM=Y-S*C0S ( A) 

RL=SQRT (ZP**2*YP**2) 

RR=SQRT (ZM**2»YM**2) 

ZPOYP=ZP/YP 
ZMOYH=ZM/YM 
PHILTLJ=ATAN (ZPOYP) 

PHIRTLJ=ATAN(ZM0YM) 

plmpi=ph1ltlj-r 

PRMPI=PrtIRTLJ-R 

COSPLI=COS(PLMPI) 

COSPRI=COS (PRMPI ) 

WNK=IS*COSPLI/RL-JS*COSPRI/RR 

RETURN 

END 


OGS 1 
OGS 2 
DGS 3 
DGS 

DGS 5 
OGS 6 
OGS 7 
OGS 8 
OGS V 
OGS 1 0 
OGS 1 1 
DGS 12 
OGS 1 3 
DGS 

DGS 15 
DGS 16 
OGS 1 7 
DGS 18 
OGS 19- 
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“DECK VLMCCIRl 

OVERLAY (miNGTL* 1*0) DGl 

program CIRCULI OGl 

OIMEMSION A0(2). 30(2)* Al(2)« 31(2)* CK2). Dl(2)* 1SUM(2). ISUMPOGl 
1(2)* ISU*<P2(2)* PPP(IOO)* WN(2)» YY<2), ZZH(50)» ZHH(IOO)* YB(50).OG1 
2 Y(IOO)* PPHK50)* XTT(50)« XTA(IOO)* CHO(IOO)* A(B*8)* C0RAG(8)* OGl 
3IP1V0T(8)* GAM(100*6)* NMA{2)« YQ(IOO). YOQ(50), YC(IOO)* YA(IOO) DGl 

COMMON /ALL/ BOT *M . BET A . &TEST « QTEST ♦ TBlSCW ( 50 ) . 0 ( AOO ) * PN ( <*00 ) *pV( 4DG1 
100) tS (400) *PS I (400) .PHI (50) ♦ZH(50) *NSSW OGl 

COMMON /ONETHRE/ T»(I$T (2) .CREE . SREF » CAVE *CL0ES ♦ STRUE . Ari ♦ ARTHUE . RTCOGl 
lOHT (2) .CONFIG.NSSWSV (2) *MSV(2) *KBOT*PLAN* IOLAN.MACH.SSVWA (50 ) .XCFWOGl 


1 

2 

3 

4 

5 

6 

7 

8 
9 


2*XCFT*YR£G( 1 *2) 

OGl 

10 

COMMON /TOTHRE/ CIR(400) 

OGl 

11 

COMMON /CCRRDO/ CHORO ( 50 ) *XTE (50 ) *KBI T.TSPAN.TSPANA 

OGl 

12 

REAL ISIGn.JSIGN 

OGl 

13 


DGl 

14 


DGl 

15 

TOLC=(80T*15.E-05) **2 

OGl 

16 


OGl 

17 


OGl 

18 

NMA ( 1 ) =NMA (2) =0 

OGl 

19 

PI=4.*ATAN ( 1 . ) 

OGl 

20 

RAO=180./PI 

OGl 

21 

B0TL=ABS(TSPAN) 

OGl 

22 

B0L=ABS(TSPANA) 

OGl 

23 

SNN=BOTL/ (2.“NSSWSV (K80T) ) 

DGl 

24 

0ELTYB=2.“SNN 

OGl 

25 

NMA (KBOT ) =B0TL/DELTY8 

DGl 

26 

NMA(KBIT)=BOl/OELTYB 

OGl 

27 

NMAX=NMA(1)*NMA(2) 

OGl 

28 

LM=1 

OGl 

29 

IF (IPLAN.EQ.2) LM=b 

OGl 

30 

IL=LM»1 

OGl 

31 

jM=LM*2 

OGl 

32 

IF (LM.EO.l) JM=IL 

OGl 

33 

IM=LM»2 

OGl 

34 

DO 1 1 = 1 *1M 

DGl 

35 

CORAG(1)=0. 

OGl 

36 

DO 1 J=1*IM 

DGl 

37 

A(I,J)=0. 

DGl 

38 

DO 2 Irl.NMAX 

OGl 

39 

DO 2 J=1*LM 

DGl 

40 

GAM(I*J)=0. 

OGl 

41 

CONTINUE 

DGl 

42 

CDRAG( ID =CLDES 

OGl 

43 

CORAG( IM) =0.0 

OGl 

44 


OGl 

45 

SCWMIN=20. 

OGl 

46 

DO 3 I=1*NSSW 

DGl 

47 

SCmMIN=AMIN1 (SCWMIN.TBLSC*MD) 

OGl 

48 

NSCWM1N=SCWMIN 

OGl 

49 

11 = 1 

OGl 

50 

DO 15 I=1.IPLAN 

OGl 

51 

B0TT=80TL 

OGl 

52 

IF (I.EO.KBIT) B0TT=80L 

OGl 

53 

IB=NSSWSV ( I ) 

OGl 

54 

1C=MSV(1)*(I-1)“MSV(2) 

OGl 

55 

IO=IC*l 

OGl 

56 

IZ=NSSWSW ( 1 ) ♦ ( l-l ) *NSSWSV(2) 

OGl 

57 

D=XCFW 

DGl 

58 

IF (I.E0.2) 0=XCFT 

OGl 

59 

AI=NSCWMIN“0*0.75 

DGl 

60 
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IMAX=INT (AI ) 

OGl 

61 


IF (O.EO.l.) 60 TO 4 

DGl 

62 


B0(I)=-1./(NSCWMIN*(1.-D)) 

DGl 

63 


AOII ) = IMAX-B0 ( I )* {'^SCWMIN*0.75)*(NSCWMIN-IMAX) 

OGl 

64 


60 TO 5 

OGl 

65 

4 

BOd )=0. 

OGl 

66 


AO(I)=IMAX 

DGl 

67 

5 

ISUM < I ) = ISUMP ( I ) *1 SUMP2 ( I ) =0 

DGl 

68 


IF (IMAX.EQ.O) GO TO 7 

DGl 

69 


DO 6 INI-1«IHAX 

OGl 

70 

b 

ISUM(I)=ISUM(I)*IN 

DGl 

71 

7 

ihm=imax*i 

DGl 

72 


IF ( IMM.6T.NSCWMIN) GO TO 9 

DGl 

73 


00 8 INsIMM.NSCWMn 

DGl 

74 


ISUMP(I)=ISUHP(I)*IN 

DGl 

75 

B 

ISUMP2 ( I ) =I SUHP2 ( I ) * IN**2 

DGl 

76 

9 

IA>4M3NHA ( I ) 

DGl 

77 


IUZ=NSSWSV(I) 

DGl 

78 


YCAT = YRF6(1 . I 1 

DGl 

79 


DO 11 J=1«IUZ 

DGl 

80 


JJ=J*(I-1)*NSSWSV(1» 

DGl 

81 


ZZH< J)=ZH(JJ) 

DGl 

82 


PPHI (J) =PhI ( JJ) 

DGl 

83 


XTT< J)=XTE(JJ) 

OGl 

84 


CIR( J)=CHORO(JJ) 

OGl 

85 


V0Q(J)=Q(I1) 

DGl 

86 


YA(JJ)=YQQ( J) 

DGl 

S7 


11=1 I*TBLSCW(JJ) 

DGl 

88 


IE=IB-J*1 

DGl 

89 


ITL=TBLSCW( IZ) 

OGl 

90 


10=IO-ITL 

OGl 

91 


IA=ID»ITL 

OGl 

92 


IF (IA.GT.IC) YCAT=YCAT-S( 10) 

OGl 

93 


IF ( IA.GT.IC) GO TO 10 

DGl 

94 


YCAT=YCAT-S( ID)-S(IA) 

OGl 

95 

10 

1Z=I2-1 

OGl 

96 


YB(IE)=YCAT 

OGl 

97 

11 

CONTINUE 

DGl 

98 


DO 12 J=1«IUZ 

DGl 

99 


JJ=J*(I-1)*NSSWSV(1) 

DGl 

100 


YC(JJ)=YB(J) 

OGl 

101 

12 

CONTINUE 

OGl 

102 


YOB = -NNA 1 1 ) *2.*SNN-SNN*YRE6( 1 1 1 ) 

DGl 

103 


00 14 K=1.IAMM 

OGl 

104 


KK=K* ( I-l )*NMA ( 1 ) 

OGl 

105 


yob=yob*oeltyb 

DGl 

106 


Y(KK)=Y0B 

DGl 

107 


CALL FTLOP (Y0B»YQ(KK) .♦1«IUZ*YB.Y0Q) 

OGl 

108 


call FTLJP (Y0R*XTA(KK) ♦♦1»IUZ»YB,XTT) 

DGl 

109 


CALL FTLOP (YOH.CHD(KK) ♦ IUZ.YB«CIR) 

OGl 

1 10 


CALL FTLOP ( YOB ♦ PP= ( KK ) * ♦ 1 ♦ I UZ » YB. PPM I ) 

DGl 

111 


call FTLOP (Y0B.ZHH<KK).*1«IUZ.Y9.ZZH) 

DGl 

112 


Bl (I )=-CM0(KK)/NSC«<MIN 

DGl 

113 


A1{I)=((XTA(KK) ♦ChO(KK) ) -0.75»B1 ( 1 ) ) *A0 ( 1 ) 

OGl 

114 


C1(I)=B0(I)«(XTA(K'<) ♦2.*CH0(KK)-1 .5*81 ( I ) ) 

OGl 

115 


D1(I)=81 (I)*BO(I) 

OGl 

116 

C 


DGl 

117 

c 

The factor b is used instead of the factor 4 to take into 

DGl 

118 

c 

ACCOUNT both sides OF THE WING 

DGl 

119 

c 


DGl 

120 
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RB=A0(I)*80(I)«ISU^P(I) 

OGl 

121 


CNNSTA=8.*SNN*C0S(ATAN(PPP<KK) ) )/SREF 

OGl 

122 


rl=cnnsta*rb 

OGl 

123 


RM=CNNSTA/CREF*(A1 ( I ) *81 ( I ) *1SUM ( I ) *01 ( I ) * ISUMP ( I ) *01 < I ) * 1 SU*^P^ ( I ) OGl 

124 


1) 

DGl 

12S 


Y8T=-YQ(KK)/B0TT 

061 

126 


SYT=SQRT ( 1 .-YBT**2) 

DGl 

127 


DO 13 JZ=1»3 

OGl 

128 


IF (IPLAN.EQ.l.ANlO.JZ.GT.l) GO TO 13 

OGl 

129 


JR=JZ*(I-1)*3 

OGl 

130 


SRL)=SYT*YBT** (2* ( JZ-1 ) ) 

DGl 

131 


GAM(KK«JR)=RB*SRU 

OGl 

132 


A(JR,IL)=A(JR«IL) ♦rl*sru 

OGl 

133 


A(JRfIM)=A(JR«IM) ♦RM*SRU 

061 

134 

13 

CONTINUE 

DGl 

135 

14 

CONTINUE 

OGl 

1 36 

IS 

CONTINUE 

OGl 

137 

c 


DGl 

138 

c 


061 

13^ 


00 16 K=1«LH 

DGl 

140 


A( IL«K ) =A (K* ID 

DGl 

141 


A( IM*K ) =A (K. IM) 

DGl 

142 

16 

CONTINUE 

061 

143 

c 


OGl 

144 

c 


OGl 

145 

c 

The -a- matrix stands for the drag matrix -cov- 

OGl 

146 

c 


OGl 

147 

c 


061 

148 


00 21 I=1,MHAX 

DGl 

149 


RPHI=ATAN(PPP(I) ) 

DGl 

130 


00 20 J=1*NMAX 

DGl 

151 


SPHI=ATAN(PPP ( J) ) 

OGl 

152 


YY ( 1 ) =YQ ( I ) -YO ( J) 

DGl 

153 


YY (2) =Y0 ( 1 ) *Y0 ( J) 

OGl 

154 


ZZ=ZHH(I)-Zhh(J) 

DGl 

155 


00 18 K=l,2 

DGl 

156 


ISIGN=JSIGN=1 . 

OGl 

157 


IF (K.E0.2) GO TO 17 

OGl 

158 


IF (YY(1) .LT.TOLC) JSIGN=-1. 

DGl 

159 


IF (YY( 1) ,LT. (-TOLC) ) ISIGN=-1. 

DGl 

IbO 

17 

YYY=YY(K) 

OGl 

161 


call DRAGSUB (RPHl.SPHI.YYY.ZZ.SNNt ISIGN. JSIGN.WN(K) ) 

061 

162 


SPHI=-SPriI 

OGl 

163 

18 

continue 

OGl 

164 


DO 19 KP=1*LM 

DGl 

165 


00 19 KG=1»LM 

DGl 

166 


A(KPtKG) =A(KP,KG) *GAM( I ,KP)*GAM( J,KG)*SNN* ( wN( 1 ) -WN(2) )/ (PI^SkEF) 

OGl 

167 

19 

CONTINUE 

DGl 

168 

20 

CONTINUE 

061 

169 

21 

CONTINUE 

OGl 

170 

C 


DGl 

1 n 



OGl 

172 


REWIND 10 

DGl 

173 


WRITE (10) ( ( A( I . J) . 1=1, JM) , J=1 ♦ JM) 

DGl 

174 


end file 10 

DGl 

175 


REWIND 20 

OGl 

176 


00 23 1=1, LM 

DGl 

177 


DO 22 J=1,LM 

DGl 

178 


XTA(J)=A(I,J)*A(J,I) 

OGl 

179 

22 

CONTINUE 

OGl 

180 


127 



APPENDIX E 


WRITE (?0) (XTA(IK) «IK=1,LM) 

23 COMTIMUE 

END FILE 20 
REWIND 20 
DO 2^ 1=1. LM 

READ (20) ( A ( I , J) » J=1 ,lM) 

2*. CONTINUE 

call SiNtO (A. JM.CORA6.1,OETER4,1PIVOT,8. ISCAlE) 

REWIND 10 

READ (10) ( ( A( I, J) . I=I , JM) , J=1 , JM) 

CD = 0. 

DO 25 1=1. LM 
DO 25 J=1.LM 

25 CD=CO*CDRAG( I) »A ( I , j) *CORAG ( J) <»2 . 

JK = 0 

DO 28 I=1.IPlAN 
B0TT=B0TL 

IF (I.E(3.KaiT) 80TT = fl0L 
KA= 1 ♦ ( I- 1 ) *N';SW5V ( 1 1 
KB=NSSWSV ( 1 ) ♦ ( l-l ) *NSSWSV (2) 

D=XCFW 

IF (I.E:0.2) D=XCFT 
DO 27 j = ->:a.kb 
rPT=-YA(J)/80TT 
SYT = SORT ( 1 .-Y8T**2) 

RJ = 0. 

DO 26 J2=1.3 
JR=JZ* ( I-l ) *3 

IF (IPLAN.EO.l.ANU.JZ.GT.l) 60 TO 26 
SRU=SYT*Y8T«« (?* ( JZ-l ) ) 

RJ=RJ*CDRAG ( JP) *SHJ 

26 continue 

NSCW=TBL5CW ( J) 

A1=NSCw*D*0.75 
IMAX=INT ( AI ) 

DO 27 K=1.NSCW 
JK = JK* 1 

E=1 . 

IF (K.GT.IMAX) E=(l.-(K-.75)/N5CW)/(l.-0) 

ClR(JK)=t*RJ 

27 CONTINUE 

28 CONTINUE 

WRITE (6.36) ClOES 
NR = 0 

DO 29 NV=1.NSSW 
NSCW=T8LSCw (NV) 

NP=NR*1 

NR=NR*NSCw 

PH I PR = AT An (PHI (NV) ) *RAD 

IF (NV.EQ. (NSSWSVU)*!) ) WRITE (6.37) 

DO 29 I=NP,NR 
PNPR=PN(I)*BETA 
PVPR=PV( I)*BETA 

PSIPR=ATAN(8ETA*TAN(PSI (I) ) )*RAO 

WRITE (6.38) PNPR«PVPR.O( I ) ,ZH(NV) ,S ( I ) .PSlPR.PHlPR.CIRl I) 

29 CONTINUE 
WRITE (6.34) 

write (6.35) CREF.CAvE«STRUE»SREF,BOT«AR»ARTRUE»MACH 
CLTOT=CMTOT=0. 

DO 31 I=1.NSSW 
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1^1 

061 

H2 

D61 


OGl 

184 

061 

IMS 

OGl 

1 86 

OGl 

187 

OGl 

18B 

OGl 

189 

OGl 

190 

OGl 

191 

OGl 

l9^ 

OGl 

193 

OGl 

194 

DGl 
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OGl 

196 

DGl 

197 

OGl 

1 98 

OGl 
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DGl 
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DGl 
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DGl 
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DGl 

205 

OGl 

206 

OGl 
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OGl 
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OGl 

209 

OGl 

210 

OGl 

211 

DGl 

212 
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213 

OGl 

214 
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061 

216 

OGl 

217 

OGl 

218 

OGl 

219 

061 
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OGl 

221 

DGl 
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OGl 
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OGl 

224 

DGl 
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DGl 
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OGl 

227 

DGl 

228 

OGl 

229 

061 

230 

OGl 
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DGl 
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OGl 
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OGl 
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DGl 
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DGl 
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OGl 

237 

DGl 
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OGl 

239 

DGl 
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30 


31 
C 

c 

32 

33 

34 


35 

36 


37 

38 

39 


IF (I.FO.l) 4RTTE (6,41) 

IF (I.EQ. (NSSwSVd)^!) ) (6,42) 

spanlo=o. 

00 30 I J=1 *NSCi*rfMIN 
IK=(I-1 )<^NSCrtMlN^I J 

SPANLO = S^AnL0^2,<^CIR ( IK) »C05 ( ATA\J(PHI (I ) ) ) 
CLTOT=CLTOTvB.^^S ( IK)<*CIR ( IK) /SPEF<^C05(ATAnmphI ( I) ) ) 


OGl 

OGl 

OGl 

OGl 

OGl 

OGl 

OGl 


IF 

IF 

IF 

IF 

IF 


( I .EQ.NSSW. and, IPlAiM.EO.2) 
( I .EQ.NSSW. and. IPLAN,E0.2) 


CONTINUE 


WRIT£(6»39) 


CMT0T=CMT0r^R.-'^5 ( IK) »CIR ( IK) <>PN ( I K ) »BET A<>COS ( AT AN ( PH I ( I ) ) ) / ( ShEF^COGI 
IREF) OGl 

Continue ogi 

WRITE (6»44) 0(IK),SPANL0 OGl 

IF ( I .EO.NSSWSV ( 1 ) ) CL1=CLT0T OGl 

( I .EO.NSSWSV (1 ) ) CM1=CMT0T OGl 

( I .EO.NSSWSV ( 1) ) WRITE (6t43) CLlfCMl OGl 

( I .EQ.NSSW. AND. IPLAN.EO.2) CL ?=CL TOT -Cl 1 OGl 

CM?=CMTOT-Cmi OGl 

WRITE (6,43) CL2,CM2 OGl 

OGl 
OGl 

CLOESf CLT0T,CMT0T,CD OGl 

OGl 

CONTINUE OGl 

Continue ogi 

RETURN OGl 

format (////4X,11H ref. chord, 6X,?5HC AVERAGE TRUE AREA ,2X,10Gl 

14HREFERENCE ARE A , 9K , 3H8/2 , 8X , 7HREF . AR,8X,7hTRUE AR , 4X , 1 1 Hm ACH NU^OGI 

2BER/) OGl 

Format (bFis.S) ogi 

format ( 1H1 ,///25X, IHXI IX, IHX, 1 1 X , 1 HY , 1 1 X , 1 HZ , I 2X , 1 HS , bX , 9 hC/ 4 SwEOGl 
1EP,4X,8HUIHE0RAL,3X, IOHGAMMA/U AT/?4X , 3HC/4 , 9X , 4H3C/4 ♦ 42 X , SHAnGLE , OG 1 
27X,5HANGlE, 4X,6HCLDES=.F7.4/) OGl 

FORMAT {/4SX,4SHSEC0ND PLANFQRM HORSESHOE VORTEX DESCRIPTIONS/) OGl 

format (17X,8F12.5) OGl 

FORMAT (/////ISXdlHCL DESIGN = , F 1 0 . 6 f 5X , 1 ?HCL COMPUTED^ ,F l 0 .6, Sx .OGl 
112HCM COMPUTED=,F10.6,bX,SHCD V=*F10.6) OGl 

FORMAT (/////1SX,7HCL 0ES=,F 1 0 .6f SX, 12HCL COMPUTED= , F 1 0 . 6 ♦ SX , 29HNOOG1 


?4l 
P42 
P4 3 

244 

245 
?46 
24 7 
247A 
2478 
?48 

249 

250 
2S0A 
?S1 
?S2 
2h2A 
?S3 

254 

255 

256 

257 

258 

259 

260 
261 
P62 

263 

264 

265 

266 

267 

268 

269 

270 

271 

272 



1 pitching moment constraint. 5X,5HC0 V=.E10.61 



DGl 

273 

41 

format 

{////40X.56HF IRST 3LANF0RM 

span 

L 0 

OGl 

274 


lA 0 I N 

G//60X . IHY. 1 IX,4HCL*C) 



OGl 

2/5 

42 

FORMAT 

(////40X.58HS ECONO PLANFORM 

span 

L 

OGl 

276 


10 A 0 I 

N G//60X. IHY, 1 1X.4HCL»C) 



OGl 

277 

43 

FORMAT 

(//50X.30HCl DEVELOPED ON THIS AnFORM= , F 1 0 

.6/ 


OGl 

278 


1 

Snx.30 HCM developed on this 3LAnFORM=.F10 

.6) 


OGl 

278A 

44 

FORMAT 

(S5XF10.5.3XF10.5) 



OGl 

279 


End 




DGl 

280- 
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*OECK VLMCC1R2 


0VERLAY(*/IN6TL*1»0» 


DG2 


program CIRCUL2 


0G2 

1 

DIMENSION A0(2)» B0(2)* Al <2) . Bl(2)« Cl(2)« 

01 (^) « 

ISOM (2). ISUMPDG2 

2 

1(2). ISUMP2(2). PP=(100)» WN(2). YY(2). ZZh(50). ZHH(IOO). YB(50).DG2 

3 

2 Y(IOO). PPHK50). XTT(SO). XTA(102). CHO(IOO). A(102.102). C0RAG(DG2 

4 

3102). 1PIVOT(102). NMA(2). YOdOO). Y00(50). 

YC(IOO) 

DG2 

5 

COMMON /ALL/ BOT .M . BFT A . PTEST . OTEST . TBlSCW ( 50 ) .0 ( ^00 ) . PN ( 400 ) .PV (4QG2 

6 

100) .5(400) ,PSI (400) .PHI (50) .ZH(50) .NSSrt 


DG2 

7 

COMMON /ONEThRE/ TWIST(2) .CREf .SREF.CA vE.CLOES.STRUE 

AkTWUt*WTC0G2 

8 

lOHT (2) .CONF IG.NSSMSV (2) .MSV (2) .K90T.PLAN. I PL AN. MACH. 

SSWWA(50) f XCFW0G2 

9 

2.XCFT.Y)^EG( 1 .2) 


0G2 

10 

COMMON /TOTHRF/ CIR(400) 


062 

11 

COMMON /CCRRDD/ CH0R0<50) .XTE (50) .K0IT.TSPAN 

.TSPANA 

0G2 

12 

real isign.jsign 


DG2 

13 



0G2 

14 



0G2 

15 

TOLC=(80T*15.E-05) »*2 


DG2 

16 



DG2 

17 



0G2 

18 

NMA(1)=NMA(2)=0 


062 

19 

PI=<*,*ATAn( 1.) 


0G2 

20 

HAO=ldO./PI 


DG2 

21 

B0TL=ABS(TSPAN) 


0G2 

22 

B0L=ABS ( TSPANA) 


062 

23 

SNN=BOTL/100. 


062 

24 

DELTYB=2.*SNN 


062 

25 

NMA (KBOT ) =aOTL/OELTYB 


062 

26 

NMA(KBIT)=B0L/0ELTYB 


062 

27 

NMAX=NMA( 1 ) .NMA(2) 


062 

28 

lm=nmax 


062 

29 

lL=LM.l 


062 

30 

jM=LM*2 


062 

31 

IF (LM.EG.NMA(l) ) JM=IL 


062 

32 

IM=LM»2 


062 

33 

DO 1 1 = 1. IM 


062 

^4 

CORAG( I ) =0. 


062 

35 

DO 1 J=1.IM 


062 

36 

A( I . J) =0. 


062 

37 

CORAG(IL) =CLDES 


062 

38 

CDRAG( IM) =0.0 


062 

39 



062 

40 

SCWMIN=20. 


062 

41 

00 2 I=1.NSSM 


062 

42 

SCWMIN=AMIN1 (SCWMIN.T8LSCW( I ) ) 


062 

43 

NSCWMIN=SCWMIN 


062 

^4 

11 = 1 


062 

45 

DO 13 I = 1.1 PL AN 


062 

46 

IB=NSSWSV ( I ) 


062 

47 

IC=MSV(1).(I-1)»MSV(2) 


062 

48 

I0=IC*1 


062 

<♦9 

IZ=NSSWSY ( 1 ) ♦ ( I-l ) »NSSWSV (2) 


062 

SO 

D=XCFW 


062 

51 

IF (I.EQ.2) 0=XCFT 


062 

52 

AI=NSCWMIN«0.0.75 


062 

53 

IMAX=INT ( AI ) 


062 

54 

IF (D.EO.l.) GO TO 3 


062 

55 

BO ( I ) =-l ./(NSCWMIN* ( l.-O) ) 


062 

56 

A0( I ) =IMAX-B0( I)*(NSCWMIN*0.75)*(NSCMMIN-IMAX) 

062 

57 

GO TO 4 


062 

58 

BO ( I ) =0 . 


062 

59 

AO ( I ) =IMAX 


062 

60 
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4 

ISUM(I)=ISJMP(I)=I5UMP2(I)=0 

DG2 

61 


IF (IMAX.Ea.O) GO TO 6 

DG2 

62 


DO 5 INsUIMAX 

DG2 

63 

D 

ISUM{I)=ISUM(I)*IN 

DG2 

64 

6 

IMM=IMAX*1 

DG^ 

66 


IF (IMM.GT.NSCWMIN) GO TO 8 

DG2 

66 


00 7 IN=IMM,NSCWMIM 

0G2 

67 


ISUMP(I)=ISUHP(I)*IN 

DG2 

68 

7 

ISUMP2( D = ISUMP2(I) ♦IN»*2 

DG2 

69 

3 

IAMM=NMA(I) 

DG2 

70 


IL»Z=NSSWSV ( I ) 

DG2 

71 


rCAT=YREG(l,I) 

DG2 

72 


DO 10 J=1.IUZ 

DG2 

73 


JJ=J* ( I-l ) *NSSWSV( 1 ) 

DG2 

74 


ZZH( J) =ZH( JJ) 

DG2 

75 


PPHI ( J)=PHI ( JJ) 

DG2 

76 


XTT ( J) =XTE ( JJ) 

DG2 

77 


CIR( J) =CHORO( JJ) 

DG2 

78 


YQO(J)=Q(II) 

DG2 

79 


1I=II*TBLSCW( JJ) 

DG2 

80 


IE=I8-J*1 

DG2 

81 


ITL=TBLSCW ( IZ) 

0G2 

82 


IO=ID-ITL 

DG2 

83 


IA*ID*ITL 

0G2 

84 


IF (IA.GT.IC) YCAT=YCAT-S(I0) 

DG2 

35 


IF (IA.GT.IC) go to 9 

DG2 

86 


YCAT=YCAT-S(IO)-StIA) 

DG2 

87 

9 

IZ=IZ-1 

DG2 

88 


YB( IE) =YCAT 

0G2 

39 

10 

CONTINUE 

DG2 

90 


00 11 J=UIUZ 

DG2 

91 


JJ=J*(I-1)*NSSWSV(1) 

0G2 

92 


YC(JJ)=Y8(J) 

DG2 

93 

11 

CONTINUE 

0G2 

94 


YOB=-NMA ( I ) *2 .*SNN-SNN*YPEG ( 1 . I ) 

DG2 

95 


00 12 K=1*IAMM 

[)G2 

96 


KK=K* ( I-l ) *NMA ( 1 ) 

DG2 

97 


Y08=Y08*OELTYB 

DG2 

98 


Y(KK)=Y0B 

DG2 

99 


CALL FTLOP ( YOB» YO (KK ) » ♦ 1 » lUZ t Y8. YOO) 

DG2 

100 


CALL FTLUP ( YOB* XT A (KK) . ♦ 1 f lUZ t YB, XTT ) 

0G2 

101 


call FTLUP (Y0B*CH0(KK).*1*IJZ.YB*CIR) 

DG2 

102 


CALL FTLUP ( YOB * PP® (KK ),♦ 1 . lUZ . YB. PPHI ) 

DG2 

103 


CALL FTLUP ( YOB . ZHH (KK ) . ♦ 1 . lUZ . YB, ZZH) 

DG2 

104 


81 (I)=-CHD(KK)/NSCi<MIN 

DG2 

105 


A1(I)=((XTA(KK) *CHD(KK) )-0.75*B1 (I ) )«A0 ( I ) 

0G2 

106 


Cl (1 )=B0 ( I )* (XTA(K-<) *2.*CH0(KK)-1.5*B1 ( I ) ) 

DG2 

107 


01 ( I ) =B1 ( I ) »B0 ( I ) 

DG2 

108 

c 


DG2 

109 

c 

The factor B is JSEO instead of the factor 4 TO TAKE 

IMTO DG2 

110 

c 

ACCOUNT both sides OF THE WING 

0G2 

111 

c 


0G2 

112 


CNNSTA=8.*SNN«C0S( ATAN(PPP(KK) ) ) /SREF 

0G2 

113 


A(kK* IL) =CNNSTA* ( AO < I ) *80 < I ) «ISUNP ( I ) ) 

DG2 

11^ 


A(KK* IM) =CNNSTA/CREF*(A1 ( I ) *9 1 ( I ) » I SUM ( I ) *C 1 ( I )*ISUMP( I ) 

(I)«ISUDG2 

115 


1MP2( I ) ) 

DG2 

1 16 

12 

CONTINUE 

DG2 

117 

13 

CONTINUE 

0G2 

118 

c 


DG2 

119 

c 


0G2 

120 


131 
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DO K=1«LM 

DG2 

121 


A(IL*K)=A(K»IL) 

0G2 

122 


A( IM»K) =A(K« IM) 

DG2 

123 

14 

CONTINUE 

DG2 

124 

C 


0G2 

125 

C 


DG2 

126 

c 

The -a- MATRIX STANDS FOR THE DRAG MATRIX -COV- 

DG2 

127 

c 


0G2 

128 

c 


0G2 

129 


00 17 I=ltLM 

DG2 

130 


RPHI=ATAN(PPP(I) ) 

0G2 

131 


CSR = A( I , IL)*SREF/(8.*SNN*C0S (RDHI ) ) 

0G2 

132 


DO 17 J=1*LM 

DG2 

133 


SPHI=ATAN(PPP(J) ) 

DG2 

1 ^4 


CSS=A( J» 1L)*SREF/(8.*SNN*C0S (SPHI ) ) 

0G2 

135 


YY(1)=YQ(I)-Y0(J) 

DG2 

136 


YY(2)=Y0(I)*Y0(J) 

0G2 

137 


ZZ=ZHH( I)-Zhh(J) 

DG2 

138 


00 16 K=l»2 

0G2 

139 


ISIGN=JSIGN=1. 

DG2 

140 


IF (K.E0.2) GO TO 15 

0G2 

141 


IF (YY( 1 ) .LT.TOLC) JSIGN=-1. 

0G2 

142 


IF (YY(l).LT.(-TOLCn ISIGN=-1. 

0G2 

143 

15 

YYY=YY (K) 

0G2 

1 <t4 


CALL DRAGSUb (RPHI .SPHI ,YYY«ZZtSNN* ISIGNt JSIGN,WN(K) ) 

0G2 

145 


SPHI=-SPHI 

DG2 

146 

16 

CONTINUE 

DG2 

147 


A( I , J)=SNN*CSR*CSS»(WN( I )-WN<2) ) / (PI*SREF) 

0G2 

148 

17 

CONTINUE 

DG2 

149 

C 


0G2 

150 

c 


0G2 

151 


REMIND 10 

D62 

152 


WRITE (10) ((A(IfJ)»I = l»JM) 

DG2 

153 


END FILE 10 

0G2 

1S4 


REMIND 20 

DG2 

155 


00 19 1=1,LM 

0G2 

156 


DO 16 J=1.LM 

DG2 

157 


XTA ( J)=A ( I « J) *A( J* I ) 

0G2 

158 

16 

CONTINUE 

0G2 

159 


WRITE (20) (XTA(IK) ,IK=1.LM) 

0G2 

160 

19 

CONTINUE 

0G2 

161 


END file 20 

0G2 

162 


REWIND 20 

0G2 

163 


DO 20 1=1. LM 

DG2 

164 


READ (20) (A(I,J).J=1.LM) 

0G2 

165 

20 

CONTINUE 

0G2 

166 


CALL SIMEQ (A. JM. CORAG. 1 .determ, IPIVOT, 102. ISCALE) 

0G2 

167 


WRITE (6. AS) 

0G2 

168 


WRITE (6.47) 

DG2 

169 


remind 10 

DG2 

170 


READ (10) ( ( A ( I . J) . 1 = 1 , JM) . J=1 , JM) 

0G2 

171 


CD=0. 

0G2 

172 


DO 21 1=1, LM 

0G2 

173 


DO 21 J=1,LM 

DG2 

174 

21 

CD=CD»CDRAG( I ) *A ( I . J)«CDRAG ( J) *2. 

0G2 

175 


DO 23 1=1. LM 

DG2 

176 


CRPHI=COS(ATAN(PPP(I) ) ) 

DG2 

177 


MNl 1=0. 

DG2 

178 


DO 22 J=1.LM 

DG2 

179 


SPHI = ATAN(PPP( J) ) 

DG2 

180 
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CSS = A ( J, ID «SRF-F / ( 3.*SNN*C0S (S=>HI ) > 

WNII=WNI I*CORAG( J) »A( I«J)*PI*SREF/(SNN*CSS*CRPHI ) 

22 CONTINUE 

write (6*46) Y { I ) *CDRAG( I ) »WNI I 

IF ( I .EO.NMA { 1 ) . ANO, IPLAN.EO.2) WPITE (6*4H) 

23 CONTINUE 

DO 26 I=1.IPLAN 
IU2 = NMA ( 1 ) 

DO 24 J=1*IU? 

JJ=J* ( 1-1 ) *NMA ( 1 ) 

ZZH( J)=Y ( JJ) 

ATT ( J)=CORAG( JJ) 

24 CONTINUE 
IUU=NSSWSV ( I ) 

DO 25 J=1«IUU 
JJ=J* ( I-l ) ®NSSWSV ( 1 ) 

CALL FTLJP ( YC( JJ) *PPP( JJ) • IUZ.ZZH.XTT ) 

25 CONTINUE 

26 CONTINUE 
JK = 0 

DO 28 1=1 . IPlAN 
KA=1 ♦(!-!) *NSSWSV ( 1 ) 

KB=NSSWSV ( 1 ) ♦ ( I-l ) *NSSwSV (2) 

D=XCFw 

IF (I.E0.2) D=XCFT 
DO 27 J=<A.KB 
NSCW=TBlSCW ( J) 

AI=NSCW»U*0.75 
IMAX=INT (A1 ) 

DO 27 K=1*NSCW 
JK=JK*1 
E = I. 

IF (K.GT.IMAX) E=(1.-(K-.75)/NSCW)/(1.-0) 

CIR ( JK ) =PPP ( J) «E 

27 continue 

28 CONTINUE 

WRITE (6*36) ClDES 
NR = 0 

DO 29 NV=1*NSSW 
NSCW=T8lSCW (NV) 

NP=NR*1 

NR=NR+NSCw 

PHIPR=ATAn (PHI (NV) )*RAO 

IF (NV.ECl. (NSSWSV( 1) ♦!) ) WRITE (6.37) 

DO 29 I=NP.NR 
PNPR=PN ( I ) *BETA 
PVPR=PV( I)*BETA 

PSIPR=ATAN(BETA*TAn(PSI (I) ) )»RAD 

write (6.38) PNPR»RVPR.0( I ) tZH(NV) »S( I ) .PSIPR.PHIPR.CIR ( I ) 

29 CONTINUE 
WRITE (604) 

WRITE (605) CREF.CAvE.STRUE.SREF.BOT.AR. ARTRUE.MACH 

cltot=cmtot=o. 

DO 31 I=1.NSSW 
IF (I.EO.l) WRITE (6.41) 

IF (I.EQ. (NSSWSVd)*!) ) WRITE (6.42) 

SPANLO=0. 

Do 30 IJ=1.NSCWMIN 
1K= ( I-l ) *NSCWM1N*I J 

SPANL0 = SPANLD*2.*CIR ( IK) »COS ( ATAnIPHI ( I) ) ) 


DG2 181 
DG2 182 
0G2 183 
D62 184 
DG2 185 
DG2 186 
062 187 
062 188 
062 189 
062 190 
062 191 
062 192 
062 193 
062 194 
062 195 
062 196 
062 197 
062 198 
062 199 
062 200 
062 pnl 
062 202 
062 203 
062 204 
062 205 
062 206 
062 207 
062 208 
062 209 
062 210 
062 211 
062 212 
062 213 
062 214 
062 215 
062 216 
062 217 
062 218 
062 219 
062 220 
062 221 
062 222 
062 223 
062 224 
062 225 
062 226 
062 227 
D62 228 
062 229 
062 230 
062 231 
062 232 
062 233 
062 234 
062 235 
062 236 
062 237 
062 238 
062 239 
062 240 
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30 


31 
C 

c 

32 

33 

34 


3b 

3b 


37 

3B 

39 

40 

41 

42 

43 

44 

45 


CLTOT = CLTOT>8.»S(IO<>CIR(IK)/S;?EF<>COS(ATAn|(PHI (I) ) ) DG2 

CMT0T = CMT0T^8.»S< I0»CIR ( IK) *Pnj ( IK ) *t)E T A«COS ( AT AN ( PH I (I ) ) ) / (SPEF<^COG2 
IPEF) DG2 

CONTINUE DG2 

WRITE (6»44) 0(IK)»SPANLD DG2 

IF (I.EO.NSSWSV(l) ) CL1=CLT0T DG2 

IF (I.FQ.NSSWSV(l) ) CM1=CMT0T OG2 

IF (I.EQ.NSSWSVU) ) WRITE (6»4l) CLl*C4l OG2 

IF (I.E0.N5SW,AND.IPLAN.E0.2) CL2=CL TOT-Cl 1 DG2 

IF (I.E0.NSSW.AND.IPLAN,EQ.2) CM2=CMTOT-Cm 1 DG2 

IF (I.EO*NSSW.AND,IPlAN.EQ.2) write (6,43) CL2,CM2 0G2 


CONTINUE 

WRITE(6,39) CEDES, CLTOT,CMTOT, CD 


OG2 
DG2 
DG2 
0G2 

CONTINUE DG2 

CONTINUE 0G2 

RETURN DG2 

format (////4X,11H ref. CHORD , 6X , 25HC AVERAGE TRUE AREA ,2X,iOG2 

14HREFERENCE ARE A , 9X , 3HB/2 t 8X , 7HREF . AR,8X,7HTRUE AR , 4X , 1 I HM ACh NUMDG2 
2BER/) DG2 

FORMAT (BF15.5) OG2 

FORMAT ( IHI ,///25X, IHXI IX, IHX, I IX, iHYf 1 IX, IhZ, 12X, 1HS»5X,9hC/4 SWEOG2 
1EP,4X,8HUIHE0RAL,3X, IOHGAMMA/U AT/24X , 3HC/4 , 9X ,4H3C/4 , 42X • BHANGLE , DG2 
27X,5HANGLE,4X,6HCLDE5=»F7.4/) DG2 

FORMAT (/45X,4SHSEC0NU PLANFORM HORSESHOE VORTEX DESCRIPTIONS/) DG2 
FORMAT (17X,8F12.5) DG2 

FORMAT (/////ISX, IIHCL DESIGN = , F 1 0 . 6 , 5X , 1 ?HCL C0MPUTED= , F 1 0 .6 , SX , 0G2 
112HCN COMPUTED=,F10.6,5X,5HCD V=»F10.6) 0G2 

FORMAT ( /////15X,7HCl OES=*F10 .6,5X , 12HCL COMPUTEO= , F 1 0 . 6 f SX , 29HNUDG2 


M 


planform 


I pitching moment constraint, 5X,5HCD V=,Fl0.b) 
FORMAT (////40X,56HF IRST ^LAnFoR 
lA 0 I N G//60X, IHY, UX,4HCL<^C) 
format (////40X,58HS ECONO 
10 A 0 I N G//60X, IHY, 1 1X,4HCL«C) 
format (//S0X,30HCl DEVELOPED ON THIS 
1 50X,30HCM developed ON THIS 

FORMAT ( 55XF10.5,3XF10 ,5) 

FORMAT (/////?X , 127HS P A N W I S E 
1 AND (NORMAL WASH) 

2H E D R A L ) ) //30X , 23H0 1 ST ANCE ALONG 
3HWN/ (U»C0S (PHI ) ) ) 


SPAN 


L 0 


span 


PLANF0RM=,Fip,6/ 
planform=,f l‘o,6) 


DG2 

OOP 

DG? 

DG2 

OG2 

DG2 

DG2 

DG2 

DO? 

DG2 


SCALE FACTOR 
/(U C0SINE([ 

PLAnFORM^SX, 7HFACTOPS,SX, 1->0G2 

002 


241 

241 A 
241B 

242 

243 

244 
;?44A 
?45 

246 
246A 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 
261 
262 

263 

264 
2bS 
266 

267 

268 

269 

270 

271 

272 
272A 

273 

274 

275 

276 

277 


46 

FORMAT 

(36XF10.5«10XF10.50XF10.5) 

DG? 

278 

47 

FORMAT 

( 10X«14hFIRST PLATFORM) 

DG? 

279 

48 

FORMAT 

( 10X« ISHSECONO PLANFORM) 

062 

280 


End 


DG2 

281 
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»OECK VLMCCIR3 

OVERLAY (WINGTL«1»0) 
program CIRCUL3 

DIMENSION A0(2)» B0(2)* Al(2)» B1 (2) ♦ Cl(2)* 01(2)* ISUM(2)* 
1(2)* ISUMP2(2)* PPP(IOO)* WN(2>* YY(2)* ZZH(50)* ZHH(IOO)* Y 
2 Y(IOO)* PPHKSO)* XTT(50). XTA(102)« CHD(IOO)* A(102*102). 
3102)* NMA(2)* YQ(IOO)* YOO(SO)* YC(IOO)* V(102*102) 

COMMON /all/ 8OT*M*BETA*PTEST,OTEST.TBlSCW(50) *Q(‘*00) *PN(‘*00 
100) *S <400 )* PS I (400) » PHI (50) *Zh(50) *NSS»( 

COMMON /ONEThRE/ TVIST (2) *CREF*SREF*CAv£*CLDES.STRiJE* ar* arth 
lOHT (2) *C0NFIG*NSSWSV <2) »MSV (2) *KBOT *PLAN * I PL AN *M ACM* SS»Mil A ( 50 
2*XCFT*YR£6( 1 *2) 

COMMON /TOTHRE/ CIR(400) 

COMMON /CCRROD/ CHORD (50 ) *XTE (50 ) *KBIT,TSPAN* TSPANA 
real ISIGN*JSIGN 


TOLC=(BOT*15.E-05)»*2 


NMA(1)=NMA(2)=0 

PI=4.*ATAN(1.) 

RAD=180./PI 

B0TL=ABS(TSPAN) 

80L=ABS(TSPAnA) 

SNN=BOTL/100, 

DELTYB*2.*SNN 

NMA(KB0T)=B0TL/DELTYB 

NMA(KBIT)=80L/DELTYB 

NMAX=NMA(1)*NMA(2) 

lm=nmax 

IL=LM*1 

JM=LM*2 

IF (LM.E'J.NMA ( 1 ) ) JM=IL 

IM=LM*2 

DO 1 I=1*1M 

CDRAG< I ) =0. 

DO 1 J=1*IM 
A(I*J)=0. 

CDRA6( IL) =CLDES 
CDRAG( IM) =0.0 

SCWMIN=20. 

DO 2 I=1*NSSW 

SC*(MIN = AM1N1 (5CWMIN*TBLSCW( I ) ) 

NSCWMINsSCrtMlN 
11 = 1 

DO 13 I=1*IPlAN 
I8=NSSWSY(I) 

IC=MSV(1)»(I-1)*MSV(2) 

IO=IC*l 

1Z=NSSWSV ( 1 ) ♦ ( I-l )*NSSWSV(2) 

0=XCFW 

IF (I.EQ.2) 0=XCFT 
AI=NSCNMIN*0*0.75 
IMAX=INT (Al ) 

IF (D.EO.l.) GO TO 3 
80(1) =-l ./ (NSCWMIN* < 1 ,-D) ) 

AO ( I ) =IMAX-B0 ( I ) * (NSCWMIN* 0.75) * (NSCWM IN- IMAX ) 

GO TO 4 

BO(I)=0. 

A0(I)=IMAX 


0G3 

HG3 

1 

I5UMPDG3 

2 

B(50) *063 

3 

CURAG(0G3 

4 

DG3 

5 

1) tPV (40G3 

6 

0G3 

7 

(Ut*PTC0G3 

8 

1 ) ♦xcrwoG3 

9 

0G3 

10 

DG3 

11 

0G3 

12 

0G3 

13 

DG3 

14 

0G3 

15 

0G3 

16 

DG3 

17 

0G3 

18 

0G3 

19 

0G3 

20 

0G3 

21 

DG3 

22 

0G3 

23 

DG3 

24 

0G3 

25 

0G3 

26 

f)G3 

27 

0G3 

28 

DG3 

29 

0G3 

30 

DG3 

31 

DG3 

32 

0G3 

33 

DG3 

34 

DG3 

35 

DG3 

36 

DG3 

37 

0G3 

38 

DG3 

39 

DG3 

40 

DG3 

41 

0G3 

42 

DG3 

43 

DG3 

44 

0G3 

4S 

0G3 

46 

DG3 

47 

0G3 

48 

DG3 

49 

0G3 

50 

DG3 

51 

DG3 

52 

0G3 

53 

063 

54 

DG3 

55 

DG3 

56 

0G3 

57 

DG3 

58 

DG3 

59 

0G3 

60 
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4 


o 

b 


7 

8 


9 

10 


11 


C 

c 

c 

c 


12 

13 

C 

c 


ISlJM(I)=ISUMP(I) = ISUHP2(I)=0 

OG3 

61 

IF ( IMAX.EO.O) GO TO 6 

0G3 

62 

DO 5 m=l*IMAX 

DG3 

63 

ISUM(I)=ISUM(I)*IN 

0G3 

64 

iMMsIMAX*! 

DG3 

65 

IF ( IMM.GT.NSCWMIN) go to 8 

DG3 

66 

DO 7 IN = IMH,NSCWHI'( 

DG3 

67 

ISUMP(I)=1SUMP(I)*IN 

OG3 

68 

ISUMP2 ( I ) = I SUMP2 ( I ) ♦ IN»*2 

OG3 

69 

IA«i1M = NMA ( I ) 

DG3 

70 

IUZ = NSSWSV( I ) 

DG3 

71 

YCAT = Y«EG(UI) 

DG3 

72 

00 10 J=1,IUZ 

OG3 

73 

JJ=J* ( I-l )*NSS«SV( 1 ) 

DG3 

74 

ZZH( J)=ZH(JJ) 

OG3 

75 

PPHI ( J) =PHI ( JJ) 

0G3 

76 

XTT{ J)=xTE( JJ) 

OG3 

77 

ClfM J)=CriORU(JJ) 

DG3 

78 

YQ0(J)=O(II) 

DG3 

79 

II = I I*TBLSCW ( JJ) 

OG3 

80 

I£=IB-J*1 

DG3 

81 

ITL=T8LSCw( IZ) 

DG3 

82 

ID=ID-ITL 

OG3 

83 

IA=IO*ITL 

0G3 

84 

IF (IA.GT.IC) YCAT=rcAT-S(IO) 

DG3 

85 

IF ( IA.GT.IC) go to 9 

DG3 

86 

YCAT=YCAT-S(ID)-S(IA) 

OG3 

87 

IZ=IZ-1 

DG3 

88 

YB( IE)=YCAT 

0G3 

89 

CONTINUE 

DG3 

90 

00 11 J=1«IUZ 

DG3 

91 

JJ=J* ( I-l ) »NSSWSV( 1 ) 

DG3 

92 

YC( JJ)=Y8(J) 

DG3 

93 

CONTINUE 

OG3 

94 

Y0fa = -NMA ( I ) *2.»SNN-SNN*YREG( 1 ♦ I ) 

DG3 

95 

00 12 Krl.lAMM 

DG3 

96 

KK=K* ( I-l ) *NMA ( 1 ) 

DG3 

97 

yob=yob*oeltyh 

0G3 

98 

Y(KK)=Y0o 

OG3 

99 

CALL FTLUP (Y0B»YQ(KK) .♦1*IUZ.Y0*YOO) 

OG3 

1 00 

CALL FTLUP (Y0B»XTA(KK) .♦l.IUZ.YB.XTT) 

0G3 

lol 

CALL FTLUP ( YOB . CHD < KK ) ♦ * 1 * I UZ , YS , C I « ) 

DG3 

102 

call FTLUP (YOH»PP=(KK) .*l»IUZ*Yg,PPHI) 

OG3 

103 

CALL FTLUP (YOBtZHH(KK) .♦l.IUZ«YB.ZZH) 

OG3 

104 

Bid) =-CHD(KK )/NSC»(MlN 

OG3 

105 

A1 ( I ) = ( (XTA (KK) ♦ChD(KK) ) - 0.75*Bl ( 1 ) ) "AO (I ) 

DG3 

106 

Cl ( I ) =B0 ( I ) " (XTA(K<) ♦2.*CH0(KK)-1 .5*B1 ( I ) ) 

DG3 

107 

0K1)=B1 (I)*B0(I) 

OG3 

108 

The factor 8 IS USED INSTEAD OF ThE FACTOR 4 

DG3 

109 

TO TAKt into DG3 

110 

ACCOUNT both sides OF THE WING 

DG3 

111 


DG3 

112 

CNNSTA=8.«SNN"C0S(ATAN(PPP(KK) ) ) /SREF 

OG3 

113 

A(KK ♦ IL ) =CNNSTA* (AO ( I ) *80 ( I ) "ISUMP ( I ) ) 

DG3 

119 

A(KK. IN) =CNNSTA/CREF* (A1 (I ) *01 ( I ) "ISUM ( I ) *C1 ( I )" 

ISUMP(I)^D1 (I)*ISUDG3 

115 

1MP2(D) 

OG3 

116 

CONT INUf 

OG3 

117 

CONTINUE 

DG3 

118 


OG3 

119 


0G3 

120 
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DO K = 1*LM 
A(1L*K)=A(K.1L) 
A( IM*K) =A (K* IM) 
I'* COI^TINUE 


the -A- MATRIX STAmOS FOR THE DRAG MATRIX -COV 


DO 17 I=1*LM 
RPHI=ATAM(PPP( I ) ) 

CSR=A ( I « IL ) *SRFF/ ( 8.*SNN*COS (RPHI ) ) 

DO 17 J=1*LM 
SPHl = ATAf'J(PPP( J) ) 

CSS=A < J* IL ) *SREF/ ( 3.*SNN*COS (St>HI ) ) 

YY(l)=YO(I)-YQ(J) 

YY(2)=YO(I) *YO(J) 

ZZ = ZHH(D-ZHH(J) 

DO 16 K=l*2 
ISIGN=JSI6H=1 . 

IF (K.EQ.2) GO TO 15 

IF (YY ( 1 ) .LT.TDLC) JSI6N=-1. 

IF (YY ( 1 ) ,LT. (-TOLC) ) ISI6N=-1. 

15 YYY=YY(K) 

CALL DRAGSUB ( RPHi tSPHI ♦ YYY . ZZ ♦ SMN* IS 1 6N . JS I GN« WN ( K) ) 
SPHI=-SPHI 

16 CONTINUE 

A(I,J)=SNN*CSR»CSS»(WN(1)-WN(2) )/(PI*S«EF) 

17 CONTINUE 
C 

C 

REWIND 10 

WRITE (10) ( (A( I .J) »1=1*JM) ♦ J=1*JM) 

END FILE 10 
REWIND 20 
DO 19 I=1«LM 
DO 18 J=1*LM 
XTA(J)=A(I«J)*A(J»n 

18 CONTINUE 

WRITE (20) (XTA(IK)fIK=l«LM) 

19 CONTINUE 
End FILE 20 
REWIND 20 

DO 20 I=1«LM 

READ (20) (A(I«J)*J=1.LM) 

20 CONTINUE 

CALL GIASOS (3fl02»102»JMtJM»A»l«CORAG*15*XTA*V*IRANK»AP«lERR) 
WRITE (6»3A) IRANKilERR 
WRITE (6*A6) 

WRITE (6*48) 

REWIND 10 

READ (10) ( (A( I ♦ J) *1 = 1* JM) ♦ J=1 ♦ JM) 

CD=0. 

DO 21 1=1. LM 
Do 21 Jsl.LM 

21 CD=C0*CDRAG( I )*A( I* J)*CDRA6( J)*2. 

DO 23 1=1. LM 
CRPHI=COS(ATAN(PPP(I) ) ) 

WNII=0. 

DO 22 J=1.LM 


DG3 121 
DG3 122 
DG3 123 
0G3 124 
0G3 125 
OG3 126 
OG3 127 
0G3 128 
OG3 129 
0G3 130 
DG3 131 
OG3 132 
DG3 1 33 
0G3 134 
OG3 135 
0G3 136 
0G3 137 
DG3 138 
DG3 139 
DG3 140 
DG3 14l 
0G3 142 
0G3 143 
DG3 144 
0G3 145 
0G3 146 
DG3 147 
0G3 148 
DG3 149 
0G3 150 
DG3 151 
0G3 152 
0G3 153 
DG3 154 
DG3 155 
0G3 156 
0G3 157 
0G3 158 
DG3 159 
0G3 160 
0G3 161 
0G3 162 
0G3 163 
DG3 164 
DG3 165 
DG3 166 
0G3 167 
DG3 168 
DG3 169 
DG3 170 
0G3 171 
DG3 172 
DG3 173 
DG3 174 
0G3 175 
0G3 176 
DG3 177 
0G3 178 
0G3 179 
0G3 180 
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SPH1 = ATAM(PPP( J) ) 

CSS=A( J, lL)*SWfc:F/(B.*SNN*COS(SPHl ) ) 

WNI 1=WNII*CDRAG( J)»A(I,J)*PI»SPEF/(SNN»CSS»CRPHI) 

22 CONTINUE 

WRITE (6*A7) Y(I) tCORAGd) ,WNII 
IF (I.EO.NMA(l) .ANO.IPLAN.EQ.2) WRITE (6.4R) 

23 CONTINUE 

DO 26 I=1»IPlAN 
IUZ = NMA ( I ) 

DO 2A J=1*IUZ 
JJ=J* ( I-l ) *NMA ( 1 ) 

ZZH(J)=Y(JJ) 

XTT ( J) =CDHA6( JJ) 

24 Continue 

IUU = NSSWSV ( I ) 

DO 25 J=1»IUU 
JJ=J*(I-1)*NSSWSV(1) 

call FTLUP (YC(JJ) .PPP(JJ) .♦l.IUZ.ZZH.XTT) 

CONTINUE 

26 Continue 
JK = 0 

DO 28 I=lfIPLAN 
KA=1MI-1)*NSSWSV(1) 

KB=NSSWSV ( 1 ) ♦ ( I-l ) »NSSWSY(2) 

D=XCFw 

IF (I.EQ.2) D=XCFT 
DO 27 J=XA,KB 
NSCW=TBlSCW( J) 

AI=nSCW*D*0.75 
IMAX=INT (AI ) 

DO 27 K=l.NSCW 
JK=JK*1 
£ = 1 . 

IF (K.GT.IBAX) E=(l.-(K-.75)/NSCW)/(l.-D) 

CIR (JK)=t>PP(J)*E 

27 CONTINUE 

28 CONTINUE 

WRITE (6»37) cedes 
NR = 0 

DO 29 NV=l«NSSw 
NSCW=TBLSCw(NV) 

NP=NR* 1 
NR=NR*NSCw 

PHIPR=ATAN (PHI (NV) ) *RAO 

IF (NV.EO. (NSSWSVd)*!) ) WRITE (6,38) 

DO 29 I=NP,NR 
PNPR=PN ( I ) *BETA 
PVPR=PV ( I ) *BETA 

PSIPR=ATAN (BETA*TAN(PSI (I ) ) ) *RAD 

write (6*39) PNPR»=VPR,Q(I)»ZH(NY) ,Sd)«PSIPR»PHIPR,CIR(I) 

29 CONTINUE 
WRITE (6*35) 

write (6*36) CREF, CAVE, STRUE,SREF,B0T, AH, ARTRUE, MACH 

cltot=cmtot=o. 

Do 31 I=1,NSSW 
IF (I.EO.l) WRITE (6,42) 

IF ( I ,EQ. (NSSWSV( 1 ) ♦! ) ) WRITE (6,43) 

SPANLD=0. 

DO 30 IJ=1,NSCWMIN 
IK= ( I-l ) "NSCwMIN*! J 


OG3 IHl 
DG3 182 
0G3 183 
063 184 
D63 185 
DG3 186 
D63 187 
OG3 188 
063 189 
0G3 190 
OG3 191 
DG3 192 
DG3 193 
DG3 194 
063 195 
0G3 196 
0G3 197 
0G3 198 
063 199 
D63 200 
063 201 
063 202 
063 203 
063 204 
063 205 
D63 206 
063 207 
063 208 
063 209 
063 210 
063 211 
063 212 
063 213 
063 214 
063 215 
063 216 
063 217 
063 218 
063 219 
063 220 
063 221 
063 222 
063 223 
D63 224 
0G3 225 
063 226 
063 227 
063 228 
063 229 
063 230 
063 231 
063 232 
063 233 
063 234 
063 235 
063 236 
063 237 
063 238 
063 239 
063 240 
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30 


31 

C 

c 

3d 

33 

34 

35 


36 

37 


38 

39 

40 

41 

42 

43 

44 

45 

46 


write ( 6*40) CLDESf CLTOTfCMTOT^CD 


SPANL0=SPANLD^2.*CIR(IK)*C0S(ATAM (PHI ( I) ) ) DG3 

CLT0T=ClT0T^8.»S( U) *CIR ( IK) /SREF<^C0S(ATAN(PHI ( I ) ) ) DG3 

CMT0T = CMT0T^8.*S( I0*CIR( IK) *PM( IK)*BETA*C0S(ATAN(?HI ( I ) ) )/ (SWEF«CDG3 
IREF) OG3 

CONTINUE DG3 

WRITE (6»45) Q(IK)*SPANLD DG3 

IF (I.EQ.NSSWSV(I) ) CLl=CLTOT DG3 

IF (I.EQ.NSSwSVd) ) CMl=CMTOT 0G3 

IF (I.EO*NSSWSV(l) ) WRITE (6^44) CLlfCHi DG3 

IF (I,EQ.NSSW.AN0.IPLAN.E0.2) CL2=CLTOT-Cl 1 DG3 

IF (I.EQ.NSSW.ANO.IPlAN.EQ.2) CM2=CMTOT-Cm I DG3 

IF (I.EQ.NSSW.AND.IPLAN.EQ.2) WRITE (6f44) CL2»CM2 DG3 

CONTINUE OG3 

OG3 
0G3 
063 

CONTINUE 0G3 

CONTINUE 0G3 

RETURN 0G3 

FORMAT ( lOXI lOdOXIlO) DG3 

FORMAT (////4XfllH REF. CHORD • 6X f 25HC AVERAGE TRUE AREA *2X,10G3 

14HREFERENCE AREA ♦ 9X ♦ 3H8/2 • 8X * 7HREF . AR,8X.7HTRUE AR ♦ 4X ♦ 1 1 HMACH NUM0G3 
2BER/) 0G3 

FORMAT (8F15.5) 063 

FORMAT ( IHI ,///25X,l NX llXf IHXfllXdHYfllX. 1 HZ* 12XdHSf5X*9HC/4 SwEDG 3 
1EP,4X*8HDIHEDRAL*3X* lOHGAMMA/U AT/24X * 3HC/4 * 9X * 4H3C/4 ♦ 42X * 5HANGLE • 063 
27X*5HANGL£*4X*6HCLDES=*F7.4/) DG3 

FORMAT (/45X,45HSEC0ND PLANFORn^ HORSESHOE VORTEX DESCRIPTIONS/) 0G3 
FORMAT (17X*8F12.5) 0G3 

FORMAT (///// 1 5X * 1 IHCL DESIGN = * F l 0 .6 * 5X * I ?hCL C0MPUTED= ♦ F i 0 . 6 * 5X * DG3 
112HCM COHPUTEO=*F10.6*5X*SHCD V=*F10.6) 063 

FORMAT (/////15X*7HCL DES= * F 1 0 .6 • SX ♦ 1 2HCL COMPUTEO= , F 1 0 .6 * 5X . 29HN00G3 


1 PITCHING MOMENT CONSTRA INT * 5X . 5HCD V=*F10.6) 

FORMAT (////40J(*56HF IRST ^LA>jFoRM SPAN 
lA 0 I N G//60XdHY* 1 1X*4HCL«C) 

FORMAT (////40X*5BHS ECONO PLANFORM SPAN 
10 A D I N G//60X* 1HY» 1 1X*4HCL<^C) 

FORMAT (//50X*30HCl DEVELOPED ON THIS ?L AnFORM= . F 1 0 . 6/ 

1 50X*30HCM DEVELOPED ON THIS ^LANFORMr » F 1 0 . 6 ) 

FORMAT ( 65XF10.5*3XF10.5) 

FORMAT (/////2X, 127HS PANWISE SCALE FACTC 
1 AND (NORMAL WASH)/(iJ * COSlNEI 


0G3 

0 DG3 
0G3 

L DG3 
0G3 
063 
063 
063 
S 0G3 

1 063 


?41 

?42 

242A 

P42B 

?43 

?44 

?45 

?45A 

?46 

?4 7 

?47A 

?48 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 
261 
262 

263 

264 

265 

266 

267 

268 

269 

270 

271 

272 
2 73 

274 
274A 

275 

276 

277 


2H E 0 R A L ) )//30X.23hOlSTANCE ALONG PL ANFORM t 5X * 7HF ACTORS f 5X * 1 5063 278 


3HWN/(U<*C0S(PHI ) ) ) 


063 279 


47 

FORMAT 

(36XF10.5*10XF10.5«3XF10.S) 

0G3 

280 

48 

FORMAT 

(10X,14HFIRST PLA*^F0RM) 

0G3 

281 

49 

FORMAT 

( lOXt ISdSECO-vtO PLANFORM) 

0G3 

282 


tMO 


DG3 

283 
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SUBHOUTINE 6IAS0S ( lOP ♦ MO ,ND . ^ , N. A , MOS . B ♦ I AC . 0 , V » IP ANK , APLUS * lERR ) 

c 

TO COMPUTE THE SINGULAR VALUE DECOMPOS I T I ON OF A REAL M X 
N MATRIX A:BY PERFORMING THE A=UOV (T) FACTORISATION, 

WITH OPTIONS FOR THE RANK,THt SINGULAR VALUtS, AN 
ORTHOGONAL BASIS FOR THE HOMUOENOUS SOLUTION , ANU THE 
PSfcUDO INVERSE OF A ANO A LEAST SOL'ARES SOLUTION FOR THE 
MATRIX problem AX=h. 


C PURPOSE 
C 
C 
C 

c 

c 

c 

C USE 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


CALL GIASOS ( IOP,Mn»NO,h,N, A,NOS,H, IAC,0, V, IRANK , APLUS, IF RR) 
lOP OPTION CODE 

IOP=l RANK WILL BE PF TURNED TO THE CALLING RROoRAR IN 

IRANK. THL ORDERLO SINGULAR VALUES WILL HF RETURNED IN 

IOP=? IN ADDITION TO ThE OPTIONS IN I0P=1 AN ORTHOGONAL 

BASIS FOR THE HOMOGENOUS SOLUTION WILL HE RETURNED In 
THE LAST N-iRANK COLUKNS OF THE V MATRIX. THE U 
transformation MATRIX WILL BF RETURNED IN MATRIX A. 


G I A 0 0 0 1 
GIA00U2 
GIA0003 
GIA0004 
GIA0005 
giaooog 
GIAOOOT 

giaooob 

GIA0009 
GlAOOlO 
GlAOOl 1 
GlAOOl? 
GIA0013 
GIA0014 
GIA0015 
GIA0016 
GIA0017 
UGI A0018 
GIA0019 
GIAOOPO 
GIA0021 
GIA0022 
GIA0023 
GIA0024 


C 

I0P = 3 

samf 

AS 

inp=2. IN addition 

THP 

least 

SuUARES 

SOLUTIONS 

01A002B 

c 


WILL 

PE 

RETURNED IN MATRIX 

B. 




G1 A0026 

c 









GIA0027 

c 

I0P = 4 

SAME 

AS 

IOR=?. IN ADDITION 

THE 

PSEUDO 

INVERSE 

WILL HE 

GIA0028 

c 


RETURNED 

IN APLUS. 





GIA0029 

c 









GIA0030 

c 

IOP = S 

SAME 

AS 

IOP=<». IN ADDITION 

THE 

LEAST 

SvlUAPES 

SOLUTIONS 

GIA0031 

c 


WILL 

HE 

RETURNfl) IN matrix 

B. 




GIA0032 


GIA0033 


c 

c 

MD 

INPUT 

INTEGER 

SPECIFING 

THE 

maximum 

ROW DIMENSION EUR A. 

6IA0034 

GIA0035 

c 

c 

NO 

INPUT 

INTEGER 

SPECIFING 

THF 

MAXIMUM 

ROW DIMENSION FOR V 

GIA0036 

GIA0037 

c 

c 

M 

INPUT 

INTEGER 

SPECIFING 

THE 

NUMBER 

OF ROWS IN A. 

GIA0038 

GIA0039 

c 

N 

INPUT 

INTEGER 

SPECIFING 

THE 

number 

OF COLUMNS IN A, 

G1A0040 


NOS 

B 


AN INPUT/OUTPUT TWO-DIMENSIONAL PEAL ARRAY WITH ROW DIMEN- 
SION MD AND COLUMN DIMENSION AT LEAST N, ON INPUTf A 
CONTAINS ThE INPUT MATRIX A WHICH IS DESTROYED. UN OUTPUT 
A CONTAINS THE ISOMETRIC MATRIX U EXCEPT WHEN IOP=l. 

NUMBER OF RIGHT HAND SIDES TO BE SOLVED, 

AN INPUT/OUTPUT TWO-DIMENSIONAL ARRAY(MD X NOS) USED FOR 


GIA0041 
GIA004? 
GIA0043 
G I AO 044 
GIA004S 
GIA0046 
GIA0047 
GIAOO40 
GIA0049 


I AC 


I0P=3 OR I0P=5, ON INPUT, B CONTAINS THE RIGHT HAND SIDES 
FOR THE SYSTEM OF EQUATIONS TO BE SOLVED. UN OUTPUT, B 
CONTAINS THE LEAST SQUARES SOLUTIONS FOR ThE EQUATIONS. 

B NEED NOT BE DIMENSIONED FOR OTHER OPTIONS. 

AN input integer SPECIFING the number of DECIMAL DIGITS OF 
ACCUPACY IN Ti^e ELEMENTS OF THE INPUT A MATPIX. THIS 
VALUE IS USEl; T(/ DETERMINE THE TEST FOR ZERO SINGULAR 
VALUES. THUS CETrW“IMNG RANK, 

IF lAC.GT.l-i The ZfcHD TFST wiLL COMPUTED USING THE 


GIA0050 
GIA0051 
GIA005? 
GI A00S3 
GIA0054 
GIA0055 
GIA0056 
&IA0057 
GIA0058 
GIA0059 
GIA0060 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


t-NOP^' OF A MULTIPLIED BY 


IF IAC.LT.13 


THE 7EPO TEST WILL BE 
E-NUBM OF A MULTIPLIED 


COMPUTED USING THE 
BY 10»*(-IAC). 


GIA0061 
GIA0062 
GIA0063 
6IA0064 
GIA0065 
GIA0066 
GIA0067 
GIA0068 

X N) WHICH CONTAINS THEGIA0069 
I0P=1. THE V MATRIX GIA0070 
CONTAIN AN ORTHOGONAL GIA0071 
N THE LAST N-IRANK GIA0072 

GIA0073 
GIA0074 
GIA0075 
GIA0076 
GIA0077 
GIA0078 
GIA0079 
GIA0080 
GIA0081 
GIA0082 
GIA0083 
GIA0084 
GIA0085 

K.GT.O IMPLIES KTH SINGULAR VALUE NOT FOUND AFTER 30 ITER. GIA0086 
K=-l implies That using THE GIVEN IAC(ACCUPACY REOUIRE- GIA0087 
mEMT)i THIS MATRIX IS CLOSE TO A mATkIX WHICh IS OF GIA0088 
LOWER RANK Than IPANK Anu IF THE ACCURACY IS GIA0089 

reduced THE Rank cf The matph 


A ONE DIMENSIONAL ARRAY OF SIZE N WHICH WILL CONTAIN THE 
ORDERED SINGULAR VALUES. 


V AN OUTPUT TWO DIMENSIONAL ARRAY (ND 

ORTHOGONAL V MATRIX EXCEPT WHEN 
UPON RETURN FROM THE SUBROUTINE WILL 
BASIS FOR THE HOMOGENOUS SOLUTIONS I 
COLUMNS FOR ALL OPTIONS EXCEPT 1 , 

IRANK RANK OF THE MATRIX A (OUTPUT) 

APLUS AN OUTPUT TWO DIMENSIONAL ARRAY {NO 
THE PSEUDO INVERSE UF MATRIX A, IF I 
4 OR 5 THIS ARRAY NEED NOT BE OIMENS 
PARAMETER MUST APPEAR IN THE CALLING 

lERR ERROR INDICATOR 

K = 0 IMPLIES normal RBTiiRN 


X M) WHICH CONTAINS 
OP DOES NOT equal 
lONEU BUT A dummy 
SEQUENCE. 


THE ACCURACY IS 
rx MAY ALSO HE REDUCED, 6IA0090 

GIA0091 


C 6IA0092 

GIA0093 

C 6IA0094 

GIA0095 
GIA0096 
G1A0097 
GIA0096 
GIA0099 
GIAOlOO 
61A0101 
GIA0102 
GIA0103 
GIA0104 
GIA0105 
GIA0106 
GIA0107 
GIA0108 
GIA0109 


LOGICAL WiTHij.wITHV 

DIMENSION A(MD,N). V ( M) . N ) ♦(! ( N ) «L ( ) 

DIMENSION H (MD*NOS) «APLUS (NU.M) 

T0L=1 .OE-60 
SIZF=0.0 
NPI=N+ 1 

COMPUTt THE E-NORM OF maTPIx A AS 7EPU TEST HOP SINGULAP VALUES 

SUM=0.0 
DO son 1=1, M 
DO 500 J=1,N 

500 SUM = SUM ♦ A(I,J)**2 
7TEST = SQPT(SUM) 


IF (IAC.GT.13) GO TO 505 GIAOllO 

7TEST = 7TEST*10.**(-IAC) GIAOlll 

GO TO 510 GIA0112 

505 7TtST = 7TEST * 2.0**(-4R) G1A0113 

7TEST =SQPT (SUM) *2.0**(-4f!) G1A0114 

C GIA0115 

510 IF (lOP.NE.l ) f,0 TO 515 GIA0116 

WITHU*. FALSE. GIA0117 

WiTHVs. FALSE. GIA0118 

GO TO 520 GIA0119 

615 W1THU=.TRUE. G1A0120 
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WiTHVr.TRUE. 

61A0121 

520 

CONTINUE 

SIA0122 


G « 0.0 

6IAQ123 


X = 0.0 

6IA0124 


DO 30 I c IfN 

6IA0125 

GIA0126 


householder reduction to bidiagonal form. 

G1A0127 

6IA0128 


Ed) = 6 

6IA0129 


S s 0.0 

GIA0130 


L = Id 

61A0131 

S1A0132 


ANNIHILATE THE I-TH COLUMN BELOW DIAGONAL. 

6IA0133 

GIA0134 


DO 3 J = IfM 

GIA0135 

3 

S = S ♦ A(JfI)*«2 

6IA0136 


G = 0.0 

6IA0137 


IF(S .LT. TOL) GO TO 10 

GIA0138 


G = SORT(S) 

GIA0139 


F = A ( I . i ) 

6IA0140 


IF(F .GE. 0.0) G = -G 

6IA0141 


H = F«G -S 

61A0142 


A(I.I) = F-G 

GIA0143 


IF(I .EG. N) GO TO 10 

GIA0144 


DO 9 J = L.N 

GIA0145 


S = 0,0 

6IA0146 


DO 7 K = I.M 

6IA0147 

7 

S = S +A (K. I ) <>A (K . J) 

bIA0148 


F = S/H 

6IA0149 


DO 8 K = I.M 

G1A0150 

8 

A(K.J) =A(K,J) ♦ F<^A(K.I) 

81A0151 

9 

CONTINUE 

GI A0152 

10 

Q( I ) = G 

G1A0153 


IF(I .EG. N) GO TO 20 

6IA0154 

GIA0155 


/iNNiMiLATE THE I-TH POw TO RIGHT OF SUPEP-OIAG. 

blAOlSft 

GIA0157 


S = 0.0 

GIA0158 


DO 11 J = L.N 

GIA0159 

11 

S = S ♦ A (I . J) A*2 

GIA0160 


G = 0.0 

6IA0161 


IF (5 .LT. TOL) GO TO 20 

GIA016? 


6 = SORT(S) 

GIA0183 


F = A(I.l4l) 

GIA0164 


IF (F .GE. 0.0) G = -G 

GIA0165 


H = F*G -S 

GIA0166 


Ad.i^n = F - G 

GIAO 167 


DO 15 J = L.N 

GIA0168 

15 

E ( J) = A (I . J) /H 

GI A0169 



DO 19 J = L.M 

GIA0170 


S = 0.0 

GIAO 171 


DO 16 K = L.K' 

GIA0172 

16 

S = S ♦ A (J.K) « A(I,K) 

GIA0173 


DO 17 K = L.N 

GIA0174 

17 

A ( J.K ) = A (J.K ) ♦ S»t (K) 

GIA0175 

19 

continue 

GIA0176 

20 

Y = ARS (0(1)) ♦ AH<; (E ( I ) ) 

GIA0177 


IF(Y .GT. SI7E) SIZE = Y 

GIA0178 

30 

CONTINUE 

GIA0179 


IE(.M0T. wlTHV) (lO TO 1 

GIA0180 
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ACCUMULATION OF FIGHT TRANSFORMATIONS. 

6IA0181 

6IA0162 


DO 40 II = I .N 

6IA0183 

6IA0184 


I s NPl - II 

6IA0185 


IF<I .EQ. N) GO TO 3R 

6IA0186 


IF(6 .EQ. 0.0) GO TO 37 

6IA0187 


H = A(I.I+1)*G 

GIA0188 


no 32 J = L.N 

GIA0I89 

3? 

V(J.I) = A(I,J)/H 

61A0190 


DO 36 J = L.N 

GIA0191 


S = 0.0 

6IA0192 


DO 33 K = L.N 

6IA0193 

33 

S = S ♦A(I.K)«V(K.J) 

6IA0194 


DO 34 K = L.N 

6IA0195 

34 

V<K.J) = V(K,J) ♦ S*V (K. I ) 

G1A0196 

36 

CONTINUE 

GIA0197 

37 

DO 3« J = L.N 

GIA0198 


V(I.J) = 0.0 

GIA0199 

38 

V ( J . I ) = 0.0 

6IA0200 

39 

V ( I . I ) = 1.0 

GIA0201 


G = F(I) 

GIA0202 

40 

L = I 

GIA0203 

41 

CONTINUE 

GIA0204 


IF(.NOT. WITHU) GO TO S3 

GIA0205 


ACCUMULATION 01- LEFT TR ANSF ORM AT I0N«. . 

GIA0206 

GIA0207 


DO 5? II = l.N 

GIA0208 

GIA0209 


I = NPl -II 

GIA0210 


L = I ♦ 1 

G1A0211 


G = 0(1) 

GIA021? 


IF (I .EQ. N) GO TO A i 

GIA0213 


DO 4? J = L.N 

GIA0214 

42 

A (I.J) = 0.0 

GIA021S 

43 

continue 

GIA0216 


IF(6 .EQ. O.C) GO TO AO 

GIA0217 


IF(I ,Eu. N) (.0 TO 4 7 

G1A0218 


H = A ( 1 , 1 ) *f- 

GIA0219 


00 46 J = L .N 

6IA0220 


S = O.n 

GIA0221 


00 4 4 K = L.r-- 

GIAO??? 

44 

S = S ♦ A(r..])AA(K.j) 

GIA0223 


F = S/r< 

G1A0224 


DO 4S K = l.N 

GIA022S 

45 

A(K,J) = A(N.J) + F*A(K.l) 

GIA0226 

46 

CONTINUE 

GIA0227 

47 

DO 48 J = I.M 

G1A0228 

48 

A ( J. I ) = A ( J, I ) /G 

GIA0229 



GO TO 51 

GIA0230 

49 

DO 50 J = I.M 

GIA0231 

50 

A ( J . I ) = 0.0 

GIA0232 

51 

A ( I . I ) = A ( I . 1 ) ♦ 1.0 

GIA0?33 

52 

CONTINUE 

GIA0234 

53 

CONTINUE 

GIA0235 



GIA0236 


DIAGONAL IZATION (E 1 U I .‘L FORM. 

GIA0237 



GIA0238 


DO 100 KK=1.N 

GIA0239 


K=NP1-KK 

GIA0240 


143 



oooo noo oooo ooo ooo 


APPENDIX E 



ITCNT=0 

6IA0241 


KP1=K+1 

GIA0242 

6IA0243 


TEST F SPLITTING. 

GIA0244 

6IA0245 

59 

CONTINUE 

6IA0246 


DO 60 LL=1.K 

6IA0247 


U=KP1-LL 

GIA0248 


IF ( (SIZE + AHS (t. (L) ) ) .EO.SI/F) GO TO 64 

6IA0249 


LM1=L-1 

GIA0250 


IF ( (SI2E + APS (0 (L"l ) ) ) .LO.SI7K ) GO TO 61 

GIA0251 

60 

CONTINUE 

GIA0252 

GIA0253 


CANCELLATION OF E(L) IF 1. .GT. 1. 

GIA0254 

GIA0255 

61 

C = 0.0 

GIA0256 


S=1 .0 

GIA0257 


L1=L-1 

6IA0258 


DO 63 I=L.K 

61 A0259 


F = S*F (I) 

6IA0260 


E (I)=C*E(I) 

GIA0261 


IF( (SIZE*ARS(F) ) .EL.SIZF) GO TO 64 

GIA0262 


G = Q ( I ) 

GIA0263 


0(1) =S0RT (F*F+G*G) 

GIA0264 


H = 0 ( I ) 

GIA0265 


C = G/H 

GIA0266 


S=-F/H 

GIA0267 


IF(.N0T.\nITHU) (j() TO 63 

GIA0268 


DO 6? J=1.M 

GIA0269 


Y=A ( J.Ll ) 

GIA0270 


Z=A(J,I) 

GIA0271 


A (J.Ll ) =Y*C*Z*S 

GIA0272 


A(J.I)= -Y«S*Z*>C 

GIA0273 

62 

CONTINUE 

61A0274 

6IA0275 

63 

CONTINUE 

6IA0276 

GIA0277 


TEST F CONVERC^ENCt . 

GIA0278 

GIA0279 

64 

7 = C5 (K) 

GIA0280 


IF(L.EG.K) GO TO 76 

GIA0281 


IFdTCNT .LE. 30) hO TO 6S 

GIA0282 


IFRR = NK 

GIA0283 


RETURN 

GIA0284 

65 

ITCNT = ITCNT ♦ 1 

GIA0285 

GIA0286 


SHIFT FROH LOWER 

GIA0287 

GIA0288 


X = 0 (L) 

GIA0289 


v = Q (K -1 ) 

GIA0290 

RsE (K-1 ) 

(iIA0291 

H=E (K) 

GIA0292 

F = ( (Y-2)«(Y + Z)* (G-h) » ((, + H) ) / (?.0«H«Y) 

6IA0293 

G = S0RT (F»F+ 1 .0 ) 

GIA0294 

IF(F.LT.O.O) 6=-G 

GIA0295 

F = ( (X-2) <» ( X + 7) +H« ( Y/ (F.G) -H) ) /X 

GIA0296 


GIA0297 


GIA0298 

NEXT OR THANSF0P^6T10N. 

GI A0299 


6IA0300 
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C=1.0 

6IA0301 


S=1.0 

GIA0302 


LP1=L^1 

GIA0303 


DO 73 I=LPl4K 

6IA0304 


G=E(I) 

6IA0305 


Y=Q(I) 

GIA0306 


H=S*G 

GIA0307 


6=C«G 

GIA0308 


2 = SORT(F*^^h<^H) 

GIA0309 


E(I-1)=Z 

GIA0310 


C=F/Z 

GIA0311 


S=H/Z 

GIA0312 


F = X»C^G<^S 

GIA0313 


Gs-X*S+G«C 

GIA0314 


H = Y*S 

GIA0315 


Y = Y*»C 

GIA0316 


IF ( .NOT.WITHV) GO TO 70 

GIA0317 


DO 68 J=1,N 

GIA0318 


X=V(JfI-l) 

6IA0319 


2=V(Jf I) 

GIA0320 


V ( Jf I-l ) =xoc + z«s 

GIA0321 


V ( Jf I ) =-X<^5^Z^C 

GIA0322 

68 

CONTINUE 

GIA0323 

GIA0324 

70 

Z = 5QPT (F»F4.H«h) 

GIA0325 


Q(I-1)=Z 

GIA0326 


C=F/Z 

GIA0327 


S = H/Z 

GIA032B 


F=sC»G + S<^Y 

GIA0329 


X=-S*G+C»Y 

GIA0330 


IF ( .NOT.WITHU) GO TO 73 

GIA0331 


DO 7? J=1,M 

GIA0332 


Y=A(J,I-1) 

6IA0333 


Z=A(J,I) 

GIA0334 


A(J,I-1)sY<^C^Z<>S 

GIA0335 


A ( J, 1 ) =-Y«S + Z^C 

GIA0336 

72 

CONTINUE 

GIA0337 

GIA0338 

GIA0339 

73 

E (L) = OdO 

GIA0340 


E(K)=F 

GIA0341 


Q (K) =X 

6IA0342 


GO TO 6P 

GIA0343 

GIA0344 


CONVERGEnCF . 

GIA0345 

GIA0346 

76 

continue 

GIA0347 


IF(Z.GE,0.O) GO TO lOQ 

GIA0348 


II 

1 

rsi 

GIA0349 



IF(*N0T.WITHV) go to 100 

GIA0350 


DO 76 vJ— 1»N 

GIA0351 

76 

V<vJ,K)=-V (J,K) 

GIA0352 

00 

CONTINUb 

6IA0353 

GIA0354 


lERR = 0 

bIA0365 


DO 280 II=2fN 

GIA0356 


1= II-l 

GIA0357 


K = I 

GIA0358 


P=Q(I) 

GIA0359 

GIA0360 
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00 se'n j=i I 

IF {0 (J) .LE.P) CjO to 2b0 
K=J 

P=0 ( J) 

?60 CONTINUE 

c 

IF (K.EQ.I) t>0 TO 2b0 
0(K) =0(1) 

0(1) = P 

C 

IF(IOP.EO.l) on TO ?hO 

c 

DO ?60 J=UN 
P= V ( J. I ) 

V(J,T)= V(J*K) 

V (J«K^ ) = P 
?P0 CONTINUE 
C 

no 2T0 J= 1 « P 
P = A ( J. I ) 

A(J*I)= A(J,K) 

A ( J,K ) = P 
270 CONTINUE 
C 

2b0 CONTINUE 
C 

J=N 

2<^0 IF (0 (U) .i)T.;TFST) (H' to jOO 
0 (J)=0.0 
J=J-1 
00 TO 2V0 
300 IRANK =J 

TFPP = 2TEST/0(U) 

IF (TFPP.GT. .0P2'-) l»P = -i 
C 

IF (lOP.LT. 3) PdiiMf 
IF(I0P.GT.3) on U 170 
DO 160 L=1.NU? 

DO 130 J=l. IRANK 
SlJP=0 . 0 
DO 120 1=1 

120 SUM =SUN + A ( I . J) OP ( I .L) 

130 F(J)= SUM/(J(J) 

C 

DO ISO K=1*N 
SUM =0.0 

DO 140 I = l,IRAr>:K 
140 SUM =SUM ♦ V (K , I ) ( 1 ) 

150 R(K,L)=SU« 


G1A0361 

GIA0362 

GIA0363 

6IA0364 

GIA0365 

6IA0366 

GIA0367 

GIA0366 

GIA0369 

GIA0370 

G1A0371 

GIA0372 

GIA0373 

GIA0374 

GIA0375 

GIA0376 

GIA0377 

GIA0378 

6IA0379 

t3IA0380 

GIA0381 

G1A0382 

6IA0383 

G1A0384 

GIA038S 

GI A0386 

t>IA0387 

(jIA038H 

GI A0389 

bl A0390 

OIA0391 

M A0392 

GI A0393 

fal A0394 

G1A0395 

C-.IA0396 

6IA0397 

GIA0398 

O1A0399 

GI A0400 

GI A0401 

GI A0402 

GIA0403 

GI A0404 

GI A0405 

GI A0406 

GIA0407 

Gl A0408 

GI A0409 


160 CONTINUE GIA0410 

RETURN GIA0411 

170 00 200 J=1.M GIA0412 

DO 190 1=1. N GIA0413 

SUM=0.0 GIA0414 

00 IPO K=l. IRANK GIA0415 

180 SUM =SUI^ V ( I ,K ) "A ( J.r ) / ; (r ) GIA0416 

190 APLUS(I.J)= SUM GIA0417 

200 CONTINUE GIA0418 

C GIA0419 

IF( lOR .EU.4) PETUPN GIA0420 
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DO 230 K^l^NOS 

GIA0421 


DO 220 I=lfN 

6IA0422 


SUMsO.O 

GIA0423 


no 210 J=lfM 

GI A0424 

210 

SUM=SUM+ APLUSdf J)<>e(J»K) 

6IA0425 

220 

E (I)=SUM 

GIA0426 


DO 225 l=l»N 

GIA0427 

225 

B(IfK)=E(I) 

GIA0428 

230 

CONTINUE 

GIA0429 


return 

GIA0430 


END 

GIA0431 
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“DECK VLMCZOC 

OVERLAY (*/INGTL*2«0) ZOC 

PROGRAM ZOCDFTm ZOC 

DIMENSION YY(2)* FV(2)« FW(2), OZDXI^^OO)* XXCC(20)» wpO(20) ZOC 

dimension X3C4(22)« AL0C(22«l)t T(4D* SS(4l*l). SSl(4l.l)« SS2U1Z0C 
1»D* S2(22»l)« S3(22«n» DELY(?2«1)« H(22)« PSUM(4l,l) ZOC 

COMMON /ALL/ HOT «M ,BET A ,PTEST ♦ QTEST . TBlSCw ( 50 ) «Q (400 ) «PN (hOO ) .PV ( 420C 
100)*S(400)»PSI (400) ♦PHI (50) f ZH(50) *NSSW ZOC 

COMMON /TOTHRE/ CIR(400) . ZOC 

COMMON /CCRROO/ CHORD ( 50 ) ♦ XTE ( 50 ) fKB I T ♦ TSP AN ♦ I SP AN A ZOC 

COMMON /1NSUR23/ APS I ♦ APHI ♦ XX ♦ Y YY ♦ ZZ ♦ SNN ♦ TOLC ZOC 

ZOC 
ZOC 

PART 3 - COMPUTE Z/C VERSUS X/C ZOC 

ZOC 
ZOC 

The tolerance set at this point in the program may need to be zoc 

CHANGED FOR COMPUTERS OTHER THAN THE CDC 6000 SERIES ZOC 

•yr\r 

£. V 

ZOC 

WRITE (6tl2) ZOC 

TOLC= (BOT*15,E-05) **2 ZOC 

IZZsl ZOC 

NNV=TBLSCW(IZZ) zoc 

DO 3 NV=ltM ZOC 

DZDX(NV)=0. ZOC 

IZ=1 ZOC 

NNN=TBLSCw(IZ) zoc 

DO 2 NN=1 fM ZOC 

APhI=ATAN (PHI ( IZ) ) zoc 

APSI=PS1(NN) zoc 

XX=PV (NV) -PN(NN) zoc 

YY (1 ) =Q (NV) -0 (NN) zoc 

YY (2) =Q (NV) »0 (NN) ZOC 

ZZ=ZH(IZZ)-Zh(IZ) zoc 

SNN=S(NN) zoc 

do 1 1=1^2 zoc 

YYY=YY(I) ZOC 

CALL INFSUB ( BOT ♦ F V ( I ) tF W ( I ) ) ZOC 

aphi=-aphi zoc 

apsi=-apsi zoc 

1 CONTINUE ZOC 

FVN=FW ( 1 ) *FW (2)-(FV(l ) ♦FV (2 ) )*phI ( IZZ) ZOC 

DZDX (NV) =OZOX (NV) ♦FVN*CIR(NN)/ 12.5663704 ZOC 

IF (NN.LT .NNN.0R.NN.F(3.M) GO T'J 2 ZOC 

IZ=IZ*1 ZOC 

NNN=NNN*THlSCW ( IZ) zoc 

2 CONTINUE zoc 

IF (NV.LT .NNV.OR.NV.EO.M) GO TO 3 ZOC 

IZZ=IZZ*1 ZOC 

NNV=NNV*THLSCW ( IZZ) zoc 

3 CONTINUE ZOC 

C ZOC 

c zoc 

C INTEGRATE DZ/DX TO OBTAIN Z/C VERSUS X/C AT THE VARIOUS Y LOCATIZOC 

C ZOC 

C ZOC 

LA=1 ZOC 

L8=0 ZOC 

DO 9 I=ltNSSW ZOC 

in=tblscw(I) zoc 

IF (I.EO.l) GO TO 4 ZOC 
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LA=LA*T8LSCW ( I-l ) 

70C 

61 

4 

LB=LB*TBLSCW ( I ) 

ZOC 

62 


00 5 J = I.A,LB 

ZOC 

63 


N=J-LA*1 

ZOC 

64 


WOU(M) =-UZ0X ( J) 

ZOC 

65 


XXCC(N) =(N-0.^S) /IV 

ZOC 

66 


K=IN*1 

ZOC 

67 


X3C4(K) =>=V ( J) *BETA 

ZOC 

68 

b 

ALOC(K) =-DZDX ( J) 

ZOC 

69 


Y=Q(LA) /aOT 

ZOC 

70 


write (6»10) Q(LA) .Y«CHORO( I ) 

ZOC 

71 


WRITE (6«13) 

ZOC 

72 


WRITE (6*17) ( WOU ( IJ) ♦ IJ=1 * I V) 

ZOC 

73 


WRITE (6*14) 

ZOC 

74 


WRITE (6*17) (XXCC(IJ)*IJ=1.IN) 

ZOC 

75 


WRITE (6*15) 

ZOC 

76 


write (6*16) 

ZOC 

77 


K1=IN*2 

ZOC 

7S 


K2=IN*1 

ZOC 

79 


AL0C(1)=AL0C(2) 

ZOC 

80 


ALOC (K1 ) =ALOC (K2) 

ZOC 

81 


X3C4( 1 ) =XTE ( I ) 

ZOC 

82 


X3C4 (Kl ) =XTE ( I ) ♦CHORD ( I ) 

ZOC 

83 


Dl=D2=0. 

ZOC 

84 


00 6 L=l*4l 

ZOC 

85 

b 

T(l)=XTE(I) ♦CHORO(I)*(L-1)*.025 

ZOC 

86 


IW = 0 

ZOC 

87 


CALL SPLINE (22*1*41, Kl, 1*41 , X 3C4, ALOC, T * A , SS * SS 1 , SS2 * S2 , S3 * OEL Y , 

HZOC 

88 


i*iw,oi,oa,i,psuM) 

ZOC 

89 


00 7 L=1 *40 

ZOC 

90 


K=42-L 

ZOC 

91 


J=41-L 

ZOC 

92 

7 

PSUM(K)=PSUM(J) 

ZOC 

93 


PSUM( 1 ) =0. 

ZOC 

94 


00 8 L=l*4l 

ZOC 

95 


K=42-L 

ZOC 

96 


X0C=1 .♦ (XTE ( I ) -T (K) )/CHORO( I ) 

ZOC 

97 


ZOC=PSUM (K) /CHORD! 1 ) 

ZOC 

98 


X = XOC*CHORO ( I ) 

ZOC 

99 


write (6*11) XOC*ZOC,X,PSUM(K) 

ZOC 

1 0 0 

b 

CONTINUE 

ZOC 

101 


WRITE (6*18) 

ZOC 

102 

9 

CONTINUE 

ZOC 

103 


RETURN 

ZOC 

104 

10 

FORMAT (35X,^HY=,F10.4* 1 1X,6HY/B/2=,F10.4, 1 1 X , 6HCH0R0= *F 1 0 .4// ) 

ZOC 

105 

1 1 

FORMAT (38XF9.4,3(5XF9.4) ) 

ZOC 

106 

Id 

format ( iHi , ssx,20hlocal elevation data///) 

ZOC 

107 

13 

FORMAT (41X,47HSL03ES,DZ/0X, AT SLOPE POINTS*FROM FRONT TO REAR/) 

ZOC 

108 

14 

FORMAT (42X,46hCORRESPONOING X/C LOCATIONS FROM FRONT TO REAR/) 

ZOC 

1 09 

lb 

FORMAT (////58X, 15HL0CAL ELEVATION//) 

ZOC 

110 

16 

FORMAT (43X,3HX/C, 1 1X*3HZ/C,8X,7HDELTA X,7X,7HDELTA Z/) 

ZOC 

111 

17 

FORMAT (5X,20F6.4) 

ZOC 

112 

18 

format (IHI) 

ZOC 

113 


end 

ZOC 

114- 
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SUBROUTIME INFSU8 ( BOT . F V I * F W I ) INF 

COMMON /INSUB23/ PSI I ♦ APH 1 1 • XXX ♦ Y YY * ZZZ t SNN ♦ TOLRNC INF 

FC=C0S(PSII ) inf 

FS=SIN(PbII) INF 

FT=FS/FC 

FPC = C0S(APHII) 

FPS=SIN(APHI I ) I*^F 

FPT=FPS/FPC inf 

fi=xxx4snn<>ft<*fpc inf 

f2=yyy4SNn<»fpc inf 

F 3 =ZZZ^SNN*FPS INF 

F4 = XXX-SNN<^FT«FPC inf 

F5 = YYY-SNN«^FPC inf 

F6 = ZZZ-SNN<>FP5 inf 

FFA= (XXX** 2 ^ ( YYY*F=>S) <^<^ 2 ^FPC** 2 * ( ( YYY*FT ) ( ZZZ/FC ) ** 2-2 . *XXX«YY iNF 

1 Y»FT ) - 2 ,^ZZZ*FPC* ( YY Y*FPS^XXX<>FT*FPS ) ) iNF 

FF 8 = (F 1 *FI^F?*F 2 ^F 3 *F 3) **.5 iNF 

FFC= (F 4 *F 4 ^F 5 *F 5 ^F 6 *F 6) **.5 INF 

FF 0 = F 5 *F 5 ^F 6 >*F 6 INF 

FFE=F 2 *F 2 ^F 3 ^»F 3 INF 

FFF= (F 1 *FPC*FT>F 2 *FPC^F 3 *FPS) /FF 9 - (F 4 *FPC*FT^F 5 *FPC^F 6 *FPS) /FFC INF 

INF 

INF 

The TOLE^^ANCE set at THIS POINT IN THE PROGRAM MAY NEED TO BE INF 

CHANGED FOR COMPUTERS OTHER THAN THE CDC 6000 SERIES INF 

INF 

INF 

IF (ABS(FFA) .LT. (B0T*15,E-S)**2) GO TO 1 INF 

FVONE=(XXX«^FPS-ZZZ«»FT<»FPC) *FFF/FFA INF 

FW0NE=( YYY*FT-XXX) «»FFF/FFA*FPC iNF 

GO TO 2 INF 

FVONE=FWONE=0. iNF 

INF 

IF ( ABS(FFD) .LT.TOLRNC) GO TO 3 INF 

FVTWO=F6«^ ( 1 .-FA/FFC) /FFO iNF 

FWTWO = -F5<»( 1 .-FA/FFC)/FFD INF 

GO TO 4 INF 

FVTWO=FWTwO=0. INF 

INF 

IF (ABS (FFE) .LT.TOLRNC) GO TO 5 INF 

fvthre=-f3<> ( i.-fi/ffb)/ffe inf 

E^^THRE = F2" ( I .-F1/FFB)/FFE INF 

GO TO 6 inf 

FVTHRE = Fi>^THRE = 0. inf 

INF 

FVI=FVOnE^FVTWO^FVTHRE inf 

fwi=fwone^fwtwo^fwthre inf 

RETURN inf 

END INF 
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6 

7 

B 

9 

10 

11 

12 

13 

14 

15 

16 


SUBROUTINE SPLINE (MNPTS^MNCVSf MHAXtNf NlCVSfM^Xf Y»T*PROXlNfSSfSSl »SSPL 
lS2*S2fS3tDELY»HfIWfDl*D2fKAB»PSUM) SPL 

DIMENSION TH(50)f DELH(50»1)» CT(50)» TH2(S0)t DELSQH(50)» ST2(50fSPL 
11) SPL 

dimension PSUM (MMAXtMNCVS) SPL 

DIMENSION X(MNPTS)* Y ( MNPTS f MNCVS ) t T(HMAX)t DELY ( MNPTS* MNCVS) t S2SPL 
1 (MNPTS^MNCYS) ♦ S3(MNPTS»MNCVS) f SSI (MMAXt MNCVS) f SS ( MMAX f M nCVS ) » HSPL 
2(MNPTS) • SS2(MMAX*MNCVS) f PROX I N ( MNCVS ) * OELSQY(50)* H2(50)f C<50)SPL 


3* 0(50) SPL 

DIMENSION DKNCVS)* D2(NCVS)» KAB(NCVS) SPL 

IF (IW) 9,1*9 SPL 


N1=N-1 

IW=2 

DO 8 K=1»NCVS 


SPL 

SPL 

SPL 


DO 2 I=1*N1 SPL 

H(I)=X(I^1)-X(I) SPL 

II = lO SPL 

OELY(IfK)=(Y(II,K)-Y(I*K) )/H(I) SPL 

C(I)=H(D SPL 

DO 3 I=2fNl .SPL 

H2(I)=(H(I-1)^H(I))*2. SPL 

OELSQY ( I ) = (DELY ( I*<)-OELY ( I-l .<) ) SPL 

CONTINUE SPL 

IF (KA0(K) ,EO.O) GO TO 4 SPL 

H2(1)=2.*H(1) SPL 

H2(N)=2.*H(N1 ) SPL 

0ELSQY(1)=6**(0ELY(1,K)-D1 (K) ) SPL 

OELSQY (N)= (02 (K) -DELY (N1 *K) )*6, SPL 

GO TO 5 SPL 

H2(l)=1.0 SPL 

H2(N)=1.0 SPL 

C(1)=0.0 SPL 

H(N1)=0.0 SPL 

DELSQY(1)=0,0 SPL 

OELSQY(N)=0*0 SPL 

CALL TRIMAT (H»H2f C»0ELSQY»0»N) SPL 

DO 6 1=1 »N SPL 

S2(IfK)=0(I) SPL 

H(n1)=C(N1) SPL 

DO 7 I=l»Nl SPL 

II=I^1 SPL 

S3(I*K)=(S2(II.K)-S2(IfK) )/H(I) SPL 

CONTINUE 5PL 

CONTINUE SPL 

J=0 SPL 

J=jO SPL 

1=1 SPL 

IF (T(J)-X(I)) 14*17,11 SPL 

IF (T(J)-X(N)) 13*15*14 SPL 

IF (T(J)-X(D) 16*17,13 SPL 

I=I^1 SPL 

GO TO 12 SPL 

CONTINUE SPL 

print 25* J SPL 

print 26* (X(I)*I=1*N) SPL 

PRINT 26* (Y(I,1)*I=1*N) SPL 

GO TO 19 SPL 

I=N SPL 

CONTINUE SPL 

IW=-I SPL 
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1=1-1 

SPL 

61 

17 

00 18 K=1»NCVS 

SPL 

t)2 


HT1=T(J)-X(I) 

SPL 

63 


II=I*1 

SPL 

64 


HT2=T(J)-X(I1) 

SPL 

65 


PR0D=HT1*HT2 

SPL 

66 


SS2(J«K)=S2(I.K)*HT1*S3<I.K) 

SPL 

67 


0ELSQS= (S2 ( I»K) »S2( II *K) *552(3 ♦,K) )/6. 

SPL 

68 


S5(J»K)=r(I*K)*HTl*DELY<I*K)*PR00*DEL53S 

SPL 

69 


551 ( J»K) =DELY(I»K) ♦ (HTl*Ht2)*DELSO5*PR0O»53( I«K)/6.0 

SPL 

70 

Id 

CONTINUE 

SPL 

71 

19 

CONTINUE 

SPL 

72 


IF (J.LT.M) GO TO 10 

SPL 

73 


M1=M-1 

SPL 

;4 


00 24 K=1*NCV5 

SPL 

75 


00 20 I=ltMl 

SPL 

76 


TH(I)=T(I*1)-T(I) 

SPL 

77 


11=1*1 

SPL 

78 


0ELH( I .K) = (5S(IItK)-S5(I»K) )/TH( I ) 

SPL 

79 


Cr (I ) =TH ( I ) 

SPL 

80 

20 

CONTINUE 

SPL 

81 


00 21 I=2«M1 

SPL 

82 


TH2(I)=( Th(I-1 ) *TH( I ) ) *2. 

SPL 

83 


0EL5QH(I)=(DELH(lf<)-OELH(I-l*K) )»6. 

SPL 

84 

21 

CONTINUE 

SPL 

85 


TH2(1)=TH2(M)=1. 

SPL 

86 


CT(1)=0 

SPL 

87 


Th(M1)=0 

SPL 

88 


DEL5QH( 1 )=DEL50H(M)=o. 

SPL 

89 


call trimat (TH»TH2*CT.0EL50H.0*M) 

SPL 

90 


DO 22 I=1.M 

SPL 

91 


ST2( ItK)=0( I) 

SPL 

92 

22 

CONTINUE 

SPL 

93 


TH(M1)=CT(M1) 

SPL 

94 


PROXIN(K) =0.0 

SPL 

95 


00 23 I=1»M1 

SPL 

96 


11=1*1 

SPL 

97 


PR0XIN(K)=P«0X1N(K)*.5«TH(I)»(55(I.K) *S5(I1,K) )-TH(I)**3*(ST2(1 

♦ K)5PL 

98 


1*ST2(I1.^) )/24. 

SPL 

99 


PSUM < 1 .K ) =PR0X1N(K) 

SPL 

100 

23 

CONTINUE 

SPL 

101 

24 

CONTINUE 

SPL 

102 


RETURN 

SPL 

103 

C 


SPL 

104 

25 

FORMAT (I4»24HTH ARGUMENT OUT OF RANGE) 

SPL 

105 

26 

FORMAT (IOFIO.3) 

SPL 

106 


END 

SPL 

107- 
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SU8R0UTIME TRIMAT ( A*B «C *0» T*N) 

OIMENSIO^J A(l)« B(l)« C(l)« 0(1)« T(l)» W(S0)* SV(SO)* 6(b0) 

This routine solves the trioiagonal (except two elements) 

W(l)=8(l) 

SV(1)=C(1)/B(1) 

G(l)=0(l)/w(l) 

NM1=N-1 
00 2 K=2*N 

^ I 

W(K)=B(K)-A(KM1)*5V<KM1) 

IF (K.EO.N) GO TO 1 
SV(K)=C(K)/W(K) 

G(K) = (0(<)-A(KM1 )*G(KM1 ) )/i«^ <K) 

CONTINUE 
T (N)=G(N) 

00 3 K=ltNMl 
KK=N-K 

T (KK) =G(f<K)-SV (KK) *T (KK^l ) 

CONTINUE 

return 

END 


TRI 

TRI 

TRI 

matriatri 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 

TRI 
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*OECK vlhcdummt 

program dummy 

UROUMB=0* 

STOP 

END 


OUM 1 
DUH 2 
OUM 3 
DUM 4- 
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ROOT-BENDING-MOMENT CONSTRAINT 


If the root bending moment is to be constrained instead of the pitching moment in 
equation (21), then it is also necessary to change from computing the Cm contributions 
(eq. (8)) to computing those of Cg. Thus the contribution to root bending moment ^ from 
the jth chordwise row would be 


'N,j 


'B.i = 


(yj - yj cos * (J, - 5^) <t>i 


b/2 


(FI) 


where is the horseshoe vortex dihedral angle. To reflect the change in the moment 
constraint which occurs in equation (21), it is necessary to rewrite the equation as 

Ns 

*^2 “ X ^k^B,k ■ ^B ~ ® (^2) 

k=l 


for the dihedral solution technique. If the constraint is that of an elliptic span loading at 
Cl d? equation (F2) becomes 


or 


(p2 


Ns 

®k^B,k 

k=l 


(-0.424414)Cl 

2 


Ns 

^2 = 2 X ^kCB,k - (-0.424414)Cl^^ 


(F3) 


(F4) 


with the number in parentheses being the fractional semispan distance of the loading cen- 
troid from the plane of symmetry. If the semispan employed in the preceding constraint 
IS different from that of the wing under consideration, as could occur in a wing with a 
winglet added on, then the fractional location must be ratioed appropriately. 

^The example of using the root-bending-moment constraint given in the text 
employed 0j = 0. 
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If there is an upper and a lower winglet on a wing to be represented, the upper wing- 
let should be defined with the wing as one planform, and the lower winglet should become 
a second planform. For this two-planform configuration, the pitching-moment constraint 
should be made inoperative and in its place, the user may want to incorporate the root- 
bending-moment constraint just described. 

The changes to the computer program listed in appendix E in order to implement 
the change to a root-bending-moment constraint are minor. The necessary change details 
will be given for the solution technique associated with the configuration having dihedral 
and a well-conditioned solution matrix. This is the part denoted by OVERLAY 1 CIRCUL2 
and on the cards by DG2. The changes are given in the order of their occurrence in the 
program: 

(1) Remove DG2 32 

(2) Replace DG2 39 with 

CDRAG(IM) = -0.424414*CLDES DG2 39 

(3) Replace DG2 115 with 

APP = ATAN(PPP(KK)) DG2 115 

and DG2 116 with 

A(KK,IM) = A(KK,IL)*((YQ(KK) - YREG(1,1))* DG2 116 

ICOS(APP) + (ZHH(KK) - RTCDHT(l))*SIN(APP))/BOT DG2 116A 

Similar changes could be made in OVERLAY 1 CIRCULl and OVERLAY 1 CIRCUL3 if 
needed. 
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Figure 1.- General layout of axis systems, elemental panels, and horseshoe vortices for 

a typical wing planform. 
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(a) a = 0.2. 

Figure 3.- Solution from two-dimensional analytical method (ref. 9) and solutions from 
present method for local slopes and elevations for various values of a. Moo = 0; 
CL^d = l-O- should be noted that the diamond symbol does not appear in the upper 
part of the figure since it is coincident with the square symbol.) 
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(a) Local slopes. 


Figure 4.- Effect of number and spanwise distribution of chordwise rows of horseshoe 
vortices on local slopes and elevations for trapezoidal wing at Cj^ ^ = 0.35 and 
Moo =0.40. A = 44.03°. 
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Figure 6.- 



Local elevations obtained from three theoretical methods for high-aspect- 


ratio wing at Cj^ d ~ 
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(b) Local elevations. 
Figure 7.- Continued. 
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Experiment (ref. 11) 
Method of ref. 1 
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(c) Lifting pressure distributions at 


Figure 7.- Concluded. 
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Figure 9.- Effect of location of chord loading change on optimum span loading, 

division, and Cq y with pitching- moment constraint for wing-canard configura- 


tion of figure 8. = 0.1; 


Moo = 0.30. 
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(c) ac = 0.6, = 0.8. 

Figure 9.- Concluded. 
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Figure 10.- Effect on span loading, Cl division, and of vertical 

separation of wing- canard configuration of figure 8 with pitching- 

moment constraint. = 0.1; a^. = 0.6; a^ = 0.8; Moo = 0.30. 
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(a) Wing. 

Figure 12.- Effect on incidence-angle distribution of vertical separation of wing-canard configuration 

of figure 8. = 0.1. 
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Figure 12.- Concluded. 
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Figure 13.- Local elevations for wing- canard configuration of figure 8 designed separately 

#-=■0.1; if 

b/2 b/2 


and in the presence of each other. Cj^ ^ = 0.2; Moo = 0.30; = O.I; = -0.676. 
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of wing 
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(b) Canard. 

Figure 13.- Concluded. 
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(a) Aerodynamic characteristics. 

Figure 15.- Effect of root-bending- moment constraint on aerodynamic characteristics, 
incidence-angle distribution, and local elevations of aspect-ratio-6.67 wing-winglet 
combination. Nc = 20; Ng = 17; Moo = 0.80. 
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(b) Incidence -angle distribution and local elevations. 
Figure 15.- Concluded. 
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